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ABSTRACT 


The  modelling  of  nonlinear  and  multidimensional  systems  from  input  and/or  output 
measurements  is  considered.  Tensor  concepts  are  used  to  reformulate  old  results  and  develop 
several  new  ones.    These  results  are  verified  through  non-trivial  computer  simulations. 

A  generalized  tensor  formulation  for  the  modelling  of  discrete-time  stationary  nonlinear 
systems  is  presented.  Tensor  equivalents  of  the  normal  equations  are  derived  and  several  efficient 
methods  for  their  solution  are  discussed.  Conditions  are  established  that  ensure  a  diagonal 
correlation  tensor  so  that  a  solution  can  be  obtained  directly  without  matrix  inversion. 

Using  a  tensor  formulation,  a  new  proof  of  the  Generalized  Lattice  Theory  is  obtained. 
Tensor  extensions  of  the  Levinson  and  Schur  algorithms  are  presented. 

New  two-dimensional  (2-D)  lattice  parameter  models  are  derived.  Using  the  tensor  form  of 
the  Generalized  Lattice  Theory  the  2-D  multi-point  error  order-updates  are  decomposed  into 
0(N  )  single  point  updates.  2-D  extensions  of  the  Levinson  and  Schur  algorithms  are  given.  The 
quarter  plane  lattice  is  considered  in  detail,  first  in  a  general  form,  then  in  forms  which  reduce  the 
computational  complexity  by  assuming  shift-invariance. 

Based  on  the  2-D  lattice,  a  new  nonlinear  lattice  model  is  developed.  The  model  is  capable 
of  updates  in  the  nonlinear  as  well  as  time  order. 


TABLE  OF  CONTENTS 

I.  INTRODUCTION   11 

A.  HISTORICAL  BACKGROUND  14 

B.  DISSERTATION  OVERVIEW  15 

II.  MATHEMATICAL  BACKGROUND  18 

A.  LINEAR  FUNCTIONALS   '. 18 

1.  Definition  2.1:  Linear  Functional  19 

2.  Theorem  2.1     19 

3.  Theorem  2.2   19 

B.  TENSORS  23 

1.  Definition  2.2:  Contravariant  Vector  23 

2.  Definition  2.3:  Covariant  Vector   24 

3.  Example  2.1    24 

4.  Definition  2.4: 'Contravariant  Tensor  of  Order  2  26 

5.  Definition  2.5:  Covariant  Tensor  of  Order  2   26 

6.  Definition  2.6:  Contravariant  Tensor  of  Order  p  26 

7.  Definition  2.7:  Covariant  Tensor  of  Order  p   26 

8.  Definition  2.8:  Mixed  Tensor  of  Order  p  27 

C.  NOTATIONAL  CONVENTIONS  28 

1.  Example  2.2    29 

2.  Example  2.3    30 

D.  BILINEAR  AND  MULTILINEAR  FORMS  31 

1.  Definition  2.9:  Bilinear  Functional    31 

2.  Example  2.4    : 32 

E.  TENSOR  OPERATIONS  33 

1.  Definition  2.10:  Tensor  Outer  Product    33 

2.  Example  2.5    34 

5 


3.  Definition  2.11:  Contraction   34 

4.  Example  2.6   35 

5.  Definition  2.12:  Tensor  Inner  Product    35 

6.  Example  2.7   36 

III.    NONLINEAR  SYSTEMS  THEORY  37 

A.  CONTINUOUS  NONLINEAR  SYSTEMS  MODELS  37 

B.  DISCRETE-TIME  NONLINEAR  SYSTEM  MODELS   39 

1.  Tensor  Formulation  42 

2.  Alternate  Tensor  Formulation  46 

3.  Example  3.1    47 

C.  DISCRETE  NONLINEAR  SYSTEM  IDENTIFICATION   48 

1.  Derivation  of  Normal  Equations  49 

2.  The  Least  Mean  Square  (LMS)  Algorithm  52 

3.  Recursive  Least  Squares  (RLS)  Algorithm  54 

4.  Simulation  Results  56 

a.  Direct  Solution  of  Normal  Equations   58 

b.  Simulation  Using  LMS  Algorithm  60 

D.  GENERALIZED  COORDINATE  SYSTEMS  60 

1.  Choices  of  Coordinate  Systems   64 

a.  Theorem  3.1    64 

b.  Proof   65 

2.  Recursive  Models  67 

a.    Example  3.2  68 

3.  Simulation  Results  71 

[V.    REVIEW  OF  LATTICE  FILTER  STRUCTURES    77 

A.    1-D  LINEAR  AUTOREGRESSIVE  LATTICE  FILTER   77 

1.    Levinson-Durbin  Algorithm  80 


a.  Theorem  4.1  Levinson's  Algorithm  (Durbin's  Form)  81 

b.  Proof  (by  Induction)  81 

2.    The  1-D  Lattice  Structure  84 

B.    GENERALIZED  ORDER  UPDATE  RECURSIONS  85 

1.  Definitions  and  Formulation  85 

2.  Generalized  Levinson  Algorithm  90 

a.  Theorem  4.2:  Generalized  Levinson  Recursion  (Regular  Form)  90 

b.  Proof  90 

c.  Theorem  4.3:  Generalized  Levinson  Recursion  (Normalized  Form)   94 

d.  Proof  94 

3.  Error  Order  Update  Recursions   96 

a.  Theorem  4.4    ^. 96 

b.  Proof  (Outline)   96 

4.  The  Generalized  Schur  Algorithm    97 

a.  Theorem  4.3:  Schur  Recursions  100 

b.  Proof   100 

5.  Synthesis  Model  101 

6.  Stochastic  Fourier  Series  Interpretation  103 

V.    TWO  DIMENSIONAL  LATTICE  STRUCTURES   107 

A.    GENERAL  FORM  OF  2-D  LATTICE  FILTER  108 

1.  Normalized  2-D  Levinson  Algorithm    110 

a.  Theorem  5.1    110 

b.  Proof  (Outline)    Ill 

2.  Normalized  2-D  Error  Order  Updates   Ill 

a.  Theorem  5.2    Ill 

b.  Proof  (Outline)   Ill 

3.  2-D  Form  of  Schur  Recursion   Ill 


a.  Theorem  5.3    112 

b.  Proof  (Outline)   112 

4.    2-D  Lattice  Structures  112 

B.  REDUCED  COMPLEXITY  2-D  LATTICE  FILTERS   116 

C.  SYNTHESIS  MODEL  122 

D.  SYSTOLIC  IMPLEMENTATIONS   126 

1.  SFG  Transformation  Procedure  127 

2.  Systolic  Implementation  of  2-D  Lattice  Filter  127 

3.  Additional  Remarks  128 

E     SIMULATION  RESULTS  129 

1.  Example  5.1    129 

2.  Example  5.2    129 

VI.  NONLINEAR  LATTICE  STRUCTURES   142 

A.  GENERAL  NONLINEAR  LATTICE  MODEL  143 

1.  Normalized  Order  Update  Recursions  146 

a.  Theorem  6.1:  Normalized  Nonlinear  Levinson  Algorithm  146 

b.  Theorem  6  2:  Normalized  Error  Order  Update  Algorithm  147 

2.  Uniqueness  Of  Lattice  Parameters  148 

a.  Theorem  6.3    148 

b.  Proof  (Outline) 148 

3.  Synthesis  Models    149 

B.  SIMULATION  RESULTS  149 

VII.  CONCLUSIONS  AND  DISCUSSION     152 

A.  Sl'MMARY  OF  NEW  RESULTS  152 

B.  FUTURE  DIRECTIONS  FOR  RESEARCH  153 

LIST  OF  REFERENCES  155 

APPENDIX  A:  ALTERNATE  PROOF  OF  THEOREM  4.2  160 


A.  DEFINITIONS  AND  FORMULATION  160 

B.  ERROR  ORDER  UPDATE  RECURSIONS  162 

1.  Theorem  4.2   162 

2.  Proof  162 

APPENDIX  B:  FORTRAN  PROGRAM  LISTINGS  164 

INITIAL  DISTRIBUTION  LIST  195 


ACKNOWLEDGEMENT 

I  owe  a  great  debt  to  Dr.  S.R.  Parker.  Without  his  constant  encouragement 
and  remarkable  insight  this  work  would  not  have  been  completed.  I  take  this 
opportunity  to  acknowledge  this  debt  and  to  thank  him. 

I  would  also  like  to  thank  the  members  of  my  Doctoral  Commitee.  In 
particular,  I  wish  to  express  my  gratitude  to  Professor  E.C.  Crittenden  for  his 
friendship  and  support,  to  Professor  P.H.  Moose  for  his  friendship  and  the 
countless  hours  which  he  spent  clarifying  many  of  my  misconceptions,  and  to 
Professor  Ziomek  for  his  careful  reading  of  this  thesis. 

Doctors  Srbijanka  Turajlic  and  Bharat  Madan  have  both  added  much  to  my 
understanding  of  Electrical  Engineering.  I  would  be  remiss  if  I  did  not  mention 
their  names. 

Professor  R.D.  Strum  has  also  contributed  substantially  to  my  education. 
However,  more  importantly,  he  gave  me  his  friendship. 

I  would  also  like  to  acknowledge  Capt(N)  J.  Dean  for  his  support  of  my 
proposal  to  remain  in  Monterey  and  persue  my  Doctorate.  In  today's  modern 
Navy  a  requirement  for  Officers  trained  to  the  Doctorate  level  exists  and  I  hope 
that  this  training  will  be  made  available  to  others  in  the  future. 

Finally.  I  must  thank  my  wife,  Kim,  and  daughters,  Kaarina.  Alexis.  Teegan. 
and  Tess.  who  endured  much  while  I  studied.  Without  this  stable  home 
environment  I  believe  this  thesis  would  not  have  been  possible. 


10 


I.    INTRODUCTION 

Linear,  time(shift)-invariant  systems  have  been  exhaustively  studied  and 
their  properties  and  behaviour  are  well  known.  These  systems  form  the 
foundation  of  all  engineering  and  scientific  disciplines.  However,  they  represent 
only  an  approximation  of  reality.  This  fact,  of  course,  does  not  diminish  their 
utility.  Mathematical  models  employing  the  assumptions  of  linearity  and  shift- 
invariance  often  provide  results  sufficiently  accurate  to  be  of  practical  use.  For  a 
large  class  of  systems,  however,  these  assumptions  cannot  be  justified  and  so 
alternate  models  must  be  used.  Mathematical  models  capable  of  representing 
nonlinearities  and  methods  for  their  identification  from  system  measurements  are 
the  major  ^topics  explored  in  this  thesis.  Due  to  the  particular  treatment  of 
nonlinear  models  chosen,  multidimensional  linear  system  modelling  is  also 
investigated.    The  assumption  of  shift-invariance  will  not  be  relaxed. 

It  will  be  useful  to  formulate  a  geometric  framework  in  which  to  solve  the 
nonlinear  modelling  problem.  The  motivation  for  this  is  simple.  The  transition 
from  physical  problems  to  geometric  ones  allows  many  diverse  phenomena  to  be 
handled  with  a  common  set  of  mathematical  tools.  This  relieves  the  burden  of 
having  to  invent  new  mathematics  for  each  new  situation.  Instead,  the  well 
understood  language  of  geometry  is  used  to  tackle  many  classes  of  problems. 
Kron  [Ref.  1:  p.  197]  states  very  clearly  the  rationale  that  allows  the  real 
physical  problem  to  be  converted  to  an  equivalent  geometric  one  in  the  following 
passage. 

A  set  of  n  equations  with  n  variables  (with  time  as  a  parameter)  mav  represent 
either  the  performance  of  a  dynamical  system  with  n  degrees  of  freedom  or  the 
motion  of  a  point  along  a  curve  located  in  an  n-dimensional  hypothetical  space 
and  expressed  along  some  frame  of  reference. 

The  basic  approach  taken  in  this  dissertation,  with  respect  to  the  modelling 
of  nonlinear  systems,  is  to  represent  them  as  a  linear  combination  of  nonlinear 
functions  of  the  data.  This  allows  linear  algorithms  to  be  applied  in  the  solution 
of  the  modelling  problem.    In  the  process  of  solving  the  nonlinear  problem  several 
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new  multidimensional  lattice  structures  are  developed. 

The  inputs  and  outputs  of  the  unknown  system  will  be  treated  as  random 
signals  (not  in  general  gaussian.)  The  problem  of  transforming  random  processes 
into  geometric  quantities  has  two  solutions.  These  are  introduced  here  in  order 
to  avoid  confusion  later  and  also  in  part  to  justify  the  tensor  formulation  that  is 
used  in  the  sequel. 

A  random  process  X  may  be  defined  as  the  assignment  of  a  function 
{x(t.u;),  teT},  to  every  outcome,  w  in  a  sample  space  Q.  Of  interest  in  this 
dissertation  is  the  case  when  T  is  a  discrete  and  finite  index  set. 

One  way  of  geometrically  visualizing  this  random  process  is  to  consider  a 
Hilbert  (or  inner  product)  space,  S,  (in  general  infinite  dimensional)  of  random 
variables,  that  is  the  vectors  or  elements  of  the  space  are  random  variables. 
.Fixing  t=t0  x(t0u;)  is  a  random  variable  and  so  is  a  vector  in  S.  The  random 
process,  X.  a  time  series  of  random  variables,  is  a  series  of  vectors  in  S.  or  a  curve 
in  S  [Ref.  2:  p.  27].  The  components  of  the  vectors  comprising  X  are  indexed  by 
the  parameter  w.  The  required  inner  product  on  this  space  is  defined  in  terms  of 
the  statistical  correlation,  ie;  E{x(tj,u;)x(uw)}.  This  approach  has  proved  highly 
successful  in  many  applications  [Ref.  3]. 

An  alternate  approach  is  to  consider  that  the  random  process  X,  consists  of 
vectors  in  a  function  space.  In  this  interpretation  the  random  process  is  an 
ensemble  of  time  functions.  {x(t,^),t£T}.  indexed  by  u.  Each  of  these  time 
functions  (generally  referred  to  as  realizations)  is  a  vector  in  the  function  space. 
There  will  exist  a  large  (in  general  infinite)  number  of  such  vectors  corresponding 
to  each  possible  outcome,  wefi.  The  components  of  the  vectors  are  indexed  by 
the  parameter  t.  There  is  no  need  to  define  a  metric  on  this  space.  Any 
expectations  that  are  required  must  be  calculated  over  the  ensemble  of  vectors. 

This  second  approach  will  be  the  one  that  is  followed  throughout  this  work. 
It  will  lead  to  many  interesting  and  novel  interpretations  of  known  algorithms 
and  also  will  be  used  to  derive  several  significant  new  results. 

While  vectors  are  sufficient  to  provide  a  complete  characterization  of 
discrete-time,    one-dimensional    linear    systems,    general    nonlinear    systems    with 
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memory  require  the  use  of  higher  order  geometric  objects  to  obtain  convenient 
descriptions.  It  is  shown  in  this  dissertation  that  a  particular  class  of  these 
geometric  objects  that  extend  vector  concepts  in  a  natural  way  and  provide  an 
ability  to  deal  with  nonlinearities  in  an  organized  fashion  are  tensors. 

The  contraversial  Sapir-Whorf  hypothesis  from  linguistics  [Ref.  4]  states 
that  the  constructs  of  a  language  define  the  boundaries  of  thought. 
Mathematics  is  a  legitimate  language.  It  is  a  well  defined  set  of  rules  used  to 
communicate  ideas.  If  the  mathematics  that  is  employed  in  the  solution  of  a 
problem  is  constrained,  then  it  is  conceivable  that  certain  solutions  may  not  be 
arrived  at,  or  even  that  the  problem  may  remain  unsolved.  In  Electrical 
Engineering,  vector  calculus  and  linear  algebra  are  the  major  mathematical  tools. 
They  are  adequate  to  explain  such  diverse  phenomena  as  the  propagation  of 
electromagnetic  waves  or  the  behaviour  of  one  dimensional  linear  systems.  More 
complex  problems  have  also  been  solved  using  this  theory  by  forcing  them  to  fit. 
but  the  notation  can  become  awkward.  Tensor  analysis  is  a  convenient 
mathematical  framework  in  which  to  deal  with  nonlinear  and  multidimensional 
signal  processing  problems.  It  provides  an  algebra  for  manipulating  objects  of 
higher  dimension  than  two,  which  is  all  that  can  easily  be  handled  using  linear 
algebra.  Importantly,  tensor  algebra  furnishes  a  system  of  notation  which  is 
powerful,  yet  compact. 

The  Electrical  Engineer's  experience  is  usually  limited  to  ordinary. 
Euclidean  geometries.  Physicists,  around  the  turn  of  the  century,  began  to 
realize  that  other  more  complex  types  of  geometries  were  equally  valid  and 
important.  In  fact.  Einstein  showed  that  the  world  we  live  in  is  neither 
Euclidean  nor  is  it  simply  three  dimensional. 

The  arguments  outlined  above  provide  the  motivation  for  this  study  of  the 
utility  of  tensor  concepts  in  Electrical  Engineering,  specifically  in  the  area  of 
discrete  signal  processing.  Although  this  work  examines  but  a  fraction  of  the 
possible  applications  in  this  field,  it  proves  that  tensors  can  lead  to  useful  results 
and  that  they  warrant  further  consideration,  particularly  in  problems  involving 
spaces  of  higher  dimensions. 

13 


A.    HISTORICAL  BACKGROUND 

Tensor  analysis  has  evolved  in  this  century,  originating  with  two  Italian 
mathematicians  in  1900;  Ricci  and  Levi-Civita.  Many  of  the  early  contributions 
are  due  to  Einstein  who  required  tensor  concepts  in  the  development  of  his 
general  theory  of  relativity.  Recent  books  on  the  subject  include  Golab,  Synge 
and  Schild,  and  Young  [Ref.  5,6,7]. 

Kron  [Ref.  l]  in  the  early  lOSO's  made  use  of  tensor  concepts  in  Electrical 
Engineering.  He  appears  to  be  first  to  do  so.  His  work  was  mainly  concerned 
with  the  analysis  and  design  of  electrical  networks  and  rotating  machinery.  Since 
that  time  there  have  been  few  papers  that  deal  with  tensors  in  the  context  of 
Electrical  Engineering. 

Volterra  [Ref.  8]  laid  the  foundation  for  nonlinear  system  analysis  in  the  late 
nineteenth  century.  He  studied  functionals  or  functions  of  functions.  He 
proposed  a  series  of  increasing  order  functionals  as  an  approximation  to  any  other 
functional.  Frechet  [Ref.  9:  p.  517]  later  showed  that  this  series  was  a  complete 
representation  and  converged  uniformly.  This  series  has  since  become  known  as 
the  Volterra  series.    We  shall  study  this  series  in  detail  in  Chapter  3. 

The  first  application  of  the  Volterra  series  to  nonlinear  systems  was  done  by 
Norbert  Wiener  in  the  1930's.  Wiener '  also  made  several  other  significant 
contributions  to  nonlinear  theory,  such  as  the  introduction  of  the  Wiener  G- 
functionals  [Ref.  10].  They  posses  the  property  of  orthogonality  when  the  system 
input  is  white  Gaussian  noise.  The  two  theories  (Volterra  and  Wiener)  form  the 
basis  of  almost  all  significant  work  to  date  on  nonlinear  systems. 

One  of  the  first  practical  methods  of  system  identification  was  proposed  by 
Lee  and  Schetzen  [Ref.  11].  Their  method  takes  advantage  of  the  orthogonality  of 
the  Weiner  G-functionals  by  employing  a  cross-correlation  technique  to  identify 
system  parameters. 

The  study  of  discrete-time  nonlinear  systems  has  gained  importance  with  the 
advent  of  the  digital  computer.  It  appears  that  the  idea  of  a  discrete  Volterra 
series  first  appeared  in  the  mid  1960's  (see  for  example  [Ref.  12].)  The  use  of 
tensor  techniques  in  the  study  of  nonlinear  systems  has  received  little  attention. 
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Sandor  and  Williamson  [Ref.  13]  made  some  use  of  them  in  the  study  of 
continuous  systems.  More  recently  Parker  and  Thomas  [Ref.  14]  proposed  the 
idea  of  using  tensor  methods  for  the  analysis  of  nonlinear  discrete-time  systems. 
Their  techniques  for  system  identification  involved  the  use  of  deterministic  input 
signals  to  extract  system  parameters. 

The  Volterra  series  is  non-recursive  and  so  a  discrete  form  cannot  represent 
an  infinite  memory  system.  This  is  equivalent  to  trying  to  represent  an  infinite 
memory  linear  system  with  a  finite  length  impulse  response.  This  can  pose 
implementation  difficulties  for  systems  with  long  memories.  One  possible  solution 
is  the  use  of  a  recursive  model.  There  has  until  very  recently  been  little  written 
about  this  because  of  the  difficulty  in  analysing  system  stability.  Parker  and 
Perry  [Ref.  15]  have  proposed  a  discrete  nonlinear  ARMA  (auto-regressive 
moving-average)  model,  however,  no  stability  implications  were  considered.  Also 
Parker.  Mayoral  and  Thomas  [Ref.  16]  proposed  an  Adaptive  Kalman  Identifier 
or  RLS  (Recursive  Least  Square)  type  algorithm  for  the  identification  of  non- 
linear ARMA  systems.  Zarzycki  and  Dewilde  [Ref.  17]  and  Zarzycki  [Ref. 
18,19,20.21]  have  proposed  a  nonlinear  lattice  structure.  Again  the  stability  of  the 
resulting  models  is  not  discussed.  Some  nonlinear  systems  are  inherently 
recursive  (eg:  the  phase  locked  loop)  so  that  this  remains  an  important  area  for 
research. 

Recently,  several  books  dealing  exclusively  with  nonlinear  systems  theory 
have  been  published.  The  book  by  Schetzen  [Ref.  9]  concerns  itself  with 
continuous  systems.  It  provides  a  very  thorough  but  readable  development  of  the 
classical  concepts.  Also  of  interest  is  a  short  appendix  outlining  the  history  of 
nonlinear  systems  theory.  A  book  by  Rugh  [Ref.  22].  is  an  important  contribution 
as  it  includes  discussions  of  discrete  theory. 

B.     DISSERTATION  OVERVIEW 

Because  the  typical  reader  of  this  dissertation  will  not  have  a  background 
which  includes  tensor  calculus  it  was  felt  that  a  chapter  covering  some 
fundamental  concepts  should  be  included.    This  material  was  considered  to  be  of 
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central  importance  to  the  work  that  followed  so  it  remains  as  a  chapter  rather 
than  being  relegated  to  an  appendix.  Readers  familiar  with  tensor  concepts  may 
wish  to  skip  most  of  Chapter  2.  although,  a  cursory  look  is  recommended  to 
ensure  that  the  notation  is  clearly  understood. 

Chapter  3  begins  with  a  review  of  the  traditional  Volterra  theory  of  nonlinear 
systems.  Both  continuous  and  discrete-time  systems  are  discussed.  The 
discrete-time  tensor  equivalent  of  the  discrete  Volterra  series  is  deduced.  An 
alternate  nonlinear  tensor  formulation  is  presented  along  with  a  discussion  of  its 
relationship  to  the  Volterra  series.  This  alternate  formulation  will  be  used  in 
most  of  the  work  that  follows. 

Next,  methods  for  the  identification  of  model  parameters  are  examined.  A 
tensor  equivalent  of  the  normal  or  Weiner-Hopf  equations  is  formulated.  Several 
recursive  in  time  algorithms  are  included  as  examples  of  the  application  of 
traditional  linear  modelling  methods  to  the  nonlinear  tensor  formulation. 
Nontrivial  numerical  simulation  results  are  included. 

The  advantages  of  using  alternate  coordinate  systems  are  then  investigated. 
It  is  shown  that  by  proper  choice  of  coordinate  systems  and  input  signals  the 
identification  process  can  be  significantly  simplified.  The  nonlinear  tensor 
formulation  is  extended  to  include  recursive  type  models.  It  is  shown  that  the 
Yule-Walker  equations  have  a  tensor  counterpart  which  can  be  solved  for  the 
model  parameters.  Several  of  the  new  results  of  this  chapter  have  already  been 
published  [Ref.  23]. 

In  Chapter  4  a  review  of  modern  lattice  theory  is  presented.  Although  the 
results  themselves  are  not  new  the  approach  is  novel.  Tensor  concepts  are  used 
to  derive  the  lattice  filters  presented  by  considering  orthogonalizing  coordinate 
transformations.  Generalized  forms  of  the  Levinson  and  Schur  algorithms  are 
also  presented  and  proven.  These  important  algorithms  are  well  known  in  linear 
matrix  theory  [Ref.  3]  and  their  generalization  in  tensor  form  is  a  significant 
result. 

Chapter  5  breaks  new  ground  by  applying  the  lattice  theory  of  Chapter  4  to 
the  problem  of  modelling  two-dimensional  data  fields.  Simplifications  due  to  an 
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assumption  of  shift-invariance  are  studied.  Several  different  configurations  are 
considered.  Simulation  results  are  included  to  prove  the  validity  of  the  theory. 
Some  implementational  aspects  of  the  algorithms  are  considered.  In  particular  a 
systolic  array  is  deduced  for  one  of  the  two-dimensional  lattice  algorithms 
presented.  This  result  demonstrates  that  the  new  algortihms  are  amenable  to 
implementation  in  dedicated  VLSI  hardware. 

In  Chapter  6  a  nonlinear  lattice  is  formulated,  again  based  upon  the  theory 
presented  in  Chapter  4  and  5.  This  is  a  new  result.  The  lattice  structure 
proposed  differs  from  those  of  previous  researchers  in  that  it  is  recursive,  not  only 
in  time  order,  but  also  in  nonlinear  order.  For  example,  one  can  obtain  the 
optimal  cubic  model  from  a  knowledge  of  the  optimal  quadratic  model.  Once 
again  nontrivial  simulation  results  are  included. 

Chapter  7  is  a  summary  of  the  new  results  presented  in  this  dissertation.  It 
draws  conclusions  about  these  results  and  outlines  some  important  unanswered 
questions  as  possible  topics  of  future  research. 

Two  appendices  are  included. 

Appendix  A  contains  an  alternate  proof  of  Theorem  4.4.  This  proof  uses  the 
Hilbert  space  formulation  described  in  this  introduction.  It  is  included  for  two 
reasons;  first,  to  illustrate  this  alternate  formulation  and  second,  to  provide 
additional  insight  into  this  theorem  which  forms  the  foundation  of  Chapters  5 
and  6. 

Appendix  B  contains  listings  of  the  FORTRAN  programs  used  in  the 
simulations  presented  in  this  thesis. 
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II.    MATHEMATICAL  BACKGROUND 

This  chapter  presents  a  brief  overview  of  the  mathematical  tools  that  are 
used  in  the  dissertation.  It  begins  with  a  discussion  of  linear,  bilinear  and 
multilinear  functionals,  and  it  is  shown  that  these  can  be  represented  by  tensors. 
Some  customary  conventions  which  simplify  notation  are  introduced,  and  some 
useful  tensor  operations  are  presented  and  discussed.  This  is  meant  to  be  an 
introduction  to  the  subject  of  tensor  analysis.  Only  those  concepts  that  will  be 
used  in  the  remainder  of  the  dissertation  are  presented.  Some  proofs  are  included 
to  act  as  examples.  Many  others  are  not  presented  here,  since  the  interested 
reader  can  find  them  in  the  references  [Ref.  5,6.7].  The  discussion  assumes  a 
thorough  knowledge  of  linear  algebra. 

A.    LINEAR  FUNCTIONALS 

We  say  that.  V,  is  a  vector  space  over  a  field  of  scalars.  F.  if  the  operations 
of  scalar  multiplication  and  vector  addition  are  defined  such  that  the  axioms  of  a 
vector  space  are  satisfied  (see  for  example  [Ref.  24.25].) 

Consider  a  vector  space.  V,  over  a  field  of  scalars.  F.  The  elements  of  V  are 
called  vectors  and  will  be  denoted  by  use  of  boldface  type,  viz  T.  If  we  restrict 
otlrselves  to  spaces  of  finite  dimension.  N,  then  we  may  write  a  vector  T  as  an 
N-tuple  of  components  and  denote  the  vector  space  by  VN.  The  components  of  a 
vector  T  will  be  denoted  by  T\  where  A  =1....N.  Writing  a  vector  as  a  set  of 
components  implies  the  existence  of  a  basis.  We  will  denote  a  basis  for  VN  as  the 

set  of  vectors  A  =  {a, aN}.    Thus  an  arbitrary  vector  T  in  VN  can  be  written 

as  a  linear  combination  of  these  basis  vectors,  viz.. 

T       £TJaA  (2.1) 

.\  - 1 

In  order  to  maintain  generality  it  is  not  necessary  to  commit  to  any  specific 
vector  space  at  this  point.    Likewise  we  allow  the  basis  to  remain  arbitrary. 
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The    following    definitions    and    theorems    are    presented    essentially    without 
proof. 

1.  Definition  2.1:  Linear  Functional 

If  VN  is  a  vector  space  over  a  field  of  scalars  F,  then,  a  linear 
transformation,  H  (the  reason  for  the  boldface  will  become  apparent  shortly,  (see 
eqn  (2.6))),  of  VN  into  F,  is  known  as  a  linear  functional  or  linear  form  on  VN. 
We  can  indicate  this  transformation  by 

H(T)  =  c        where  ceF  and  TeVN  (2.2) 

2.  Theorem  2.1 

The  set  of  all  linear  functionals  on  VN  forms  a  vector  space  of  the  same 
dimension  as  VN.  This  space  is  known  as  the  dual  vector  space  and  is  indicated 
by  VN. 

3.  Theorem  2.2 

If  VN  is  a  vector  space  over  the  field  of  scalars  F,  with  basis  A  = 
{aj,  ...,  aN},  then  the  set  of  linear  functionals  A  =  {b1  ....  bN},  (defined  so  that  the 
A-th  functional.  bA.  operating  on  an  arbitrary  element  of  VN.  say  T.  yields  the  A- 
th  component  of  T.  namely  T*  ).  form  a  basis  for  VN.  The  defining  property  can 
be  expressed  mathematically  as 

bA(T)  =  TA        for  all  A  6  1 N  (2.3) 

N 
where  T  =    J]T  a^ 
A=  I 

We  call  the  functional  bA  the  A-th  coordinate  function  since  when  applied 
to  a  vector  T  it  yields  the  A-th  coordinate,  namely  T\  The  set  of  these  linear 
functionals,  b\  A  e     {1 N}  comprising  A  is  known  as  The  dual  basis  of  A. 

It  is  interesting  to  note  that  this  choice  of  basis  for  the  dual  space  leads 
to  the  property  that 

J  ]    when  A-  fj  ,<->   <  \ 

b'(aj  =  S\  =  in       .       .   ,  I2-4) 

v    M/  M         10    whenA^// 

To  show  this  we  proceed  as  follows:  from  (2.1)  it  follows  that 
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/n 
bA(T)  =  bA   £T"a„ 

V  =  1 


=    i>A(T"a«) 


(2.5a) 


(2.5b) 


Since  bA  is  a  linear  functional  (2.5b)  can  be  written  as, 


bA(T)  =    £T*bA(a„) 


(2.5c) 


However,  from  (2.3)  it  is  known  that  bA(T)  =  TA.    Thus  (2.5c)  implies  (2.4). 

The  existence  of  a  basis  for  the  dual  space  implies  that  any  vector  (linear 
functional)  H  in  VN  may  be  written  uniquely  as  a  linear  combination  of  the 
elements  of  the  dual  basis.  Therefore,  any  linear  functional  can  be  represented 
uniquely  by  an  N-tuple  of  components.    Thus 


H=   £HAb^ 


A=l 


From  (2.6)  one  can  write 


H(T)  =    £H,bA(T) 
Using  (2.3),  (2.7a)  becomes 


H(T)  =    1]HaTa 

A=l 

Alternatelv.  in  matrix  form,  if 


H  =   iH,  H- 


HN 


(2.6) 


r2.7a: 


(2.7b) 


'2.81 


T  _  [Ti  T2  .  .  .   tn  T  (2.9) 

then  (2.6)  can  be  written  as: 

H(T)  -   HT.  (2.10) 

We  notice  that  the  two  vectors.  H  and  T.  are  defined  relative  to  different 
basis   and   that   they   belong   to   different   vector   spaces.     The   vector   H.   defined 
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according  to  the  dual  basis  A  is  called  a  covariant  vector.  The  vector  T.  defined 
according  to  the  regular  basis.  A.  is  known  as  a  contravariant  vector.  As  a 
convention,  whenever  a  subscript  is  used  to  index  the  components  of  a  vector  it 
is  understood  that  the  vector  is  being  expressed  according  to  the  dual  basis,  A,  in 
the  dual  vector  space  VN,  and  is  a  covariant  vector.  Similarly,  when  a 
superscripted  index  is  used  the  components  are  assumed  to  represent  a 
contravariant  vector  in  the  vector  space  VN  according  to  the  regular  basis  A. 

Equation  (2.7b)  represents  what  we  normally  think  of  as  a  vector  inner 
product.  We  usually  do  not  think  of  the  two  vectors  as  coming  from  different 
vector  spaces.  The  reason  for  this  is  that  in  the  familiar  rectangular  cartesian 
system  of  coordinates,  the  regular  and  dual  basis  are  identical  and  so  there  is  no 
need  to  differentiate  between  covariant  and  contravariant  vectors.  In  other 
systems  of  coordinates  the  distinction  must  be  made  in  order  that  the 
relationships  have  meaning.  To  perform  a  vector  inner  product,  one  vector  must 
be  covariant  and  the  other  must  be  contravariant.  For  example  consider  the 
vector  T  as  illustrated  in  Figure  2.1.  It  has  components  l  3T  with  respect  to 
basis  {e!,e2}  and  components    -l    2}T  with  respect  to  basis  {er.e2-}. 

A  measure  of  the  length  of  vector  T  in  the  rectangular  coordinate  system. 
{e1,e2},  can  be  computed  from  the  expression 


N 

£TAT; 

A=  1 


1/  2 


(2.11a) 


1/  2 


(2.11b) 


10 


(2.11c) 


The  answer  i^  correct  because  in  the  rectangular  ("artesian  coordinate  system 
there  is  no  need  to  distinguish  between  covariant  and  contravariant  vectors,  ie: 
T^  =  TA.  However,  an  expression  similar  to  (2.11c)  in  the  oblique  coordinate 
system,  {ex   e2  }.  does  not  yield  a  measure  of  the  length  of  the  vector. 
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Figure  2.1:  Vector  T  expressed  according  to  the  two  bases  {e,.e2}  and  {er.e2-}. 


22 


N  r  ni/2 

^TATA=       (_1)2^22  (212a) 

A  =  l  L  J 


51    - 


(2.12b) 


The  answer  is  incorrect  since  both  vectors  in  expression  (2.12)  are  contravariant. 
To  obtain  a  correct  answer,  one  of  the  vectors  would  have  to  be  made  covariant. 
This  involves  introduction  of  a  metric  tensor.  The  interested  reader  can  find  a 
discussion  of  this  concept  in  [Ref.  6].  In  the  case  of  rectangular  coordinate 
systems  this  metric  tensor  is  the  identity  matrix.  Our  intent  here  was  only  to 
indicate  that  in  any  system  other  than  rectangular  Cartesian,  strict  attention 
must  be  paid  to  the  character  of  the  vectors. 

B.     TENSORS 

In  the  previous  section  the  concept  of  covariant  or  contravariant  vectors  has 
been  established.  In  this  section  a  more  formal  definition  of  these  quantities  is 
presented. 

Suppose    we    have    N    variables    x\  x2 xN.    then    a    set    of  values    of   these 

variables  is  called  a  point.  The  variables  themselves  are  called  coordinates  (or 
components.)     The  totality  of  all  points,  as  each  of  the  variables  (coordinates) 

x\  A  -  l N,  vary  over  their  entire  specified  range,  constitutes  an  N-dimensional 

space,  denoted  by  VN. 

1.      Definition  2.2:  Contravariant  Vector 

A  contravariant  vector  T.  is  defined  on  the  basis  of  the  transformation  of 
its  components  upon  transition  from  one  coordinate  system  to  another.  For 
coordinate  system  (A)  the  components  of  T  are  an  X-tuple  of  numbers  designated 
as 

TA    (A  =  ] \) 

Upon  transition  to  another  coordinate1  system  (A'),  if  the  components  of  T 
transform  according  to  the  rule 
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TA'  =   E^T^  (2.13) 

where  xA  and  xA'   define   the   coordinates  of  a   point    in  the  old   (A)    and  new   (a') 
coordinate  systems,  then  T  is  said  to  be  a  contravariant  vector. 

2.  Definition  2.3:  Covariant  Vector 

A  covariant  vector  U,  is  defined  on  the  basis  of  the  transformation  of  its 
components  upon  transition  from  one  coordinate  system  to  another.  For 
coordinate  system  (a)  the  components  of  U  are  an  N-tuple  of  numbers  designated 
as 

U,      (A  =  1,  ....  N) 

Upon    transition    to    another    coordinate    system    (A'),    if  the    components    of  U 
transform  according  to  the  rule 

Ur=E^UA  (2.14) 

X=  i  ox 

then  U  is  said  to  be  a  covariant  vector. 

g     A  ' 

In  equation  (2.13)  the  quantity  represents  the  partial  derivative  of 

the  new   (primed)  coordinates  with  respect  to  the  old  coordinates.    Similarly,  in 

ft   x 
equation  (2.14)  — —  represents  the  partial  derivative  of  the  old  coordinates  with 
3x 

respect    to    the    new.    primed,    coordinates.     In    general    these    quantities    can    be 
arranged  into  a  two-dimensional  matrix  of  numbers. 

3.  Example  2.1 

Consider  the  following  parametric  description  of  a  curve: 

x1  =  f,(t)  12. 15a) 

x2  =  f2(t)  (2.15b) 

xs  -  f,(t)  (2.15c) 

We  consider  the   three   quantities  x1.  x:.  x\  to   be   the  coordinates  of  a   point   or 
equivalently  the  components  of  a  vector  in  a  three  dimensional  space.    We  leave 


24 


the  basis  arbitrary.  Indeed  the  equations  we  write  will  be  true  regardless  of  the 
choice  of  basis.  This  is  the  inherent  advantage  of  tensor  analysis  since  it  allows 
expressions  to  be  written  which  are  invariant  with  respect  to  the  coordinate 
system. 

We  can  now  define  new,  primed,  coordinates  as  functions  of  the  old 
coordinates.  For  example,  if  we  arbitrarily  choose  the  following  coordinate 
transformation 


x*'  =  MxV.x') 


(2.16a) 


x2    =  g2(x\x2.x3) 


(2.16b) 


x"    =    g3(x  ,X  .X  ) 


—  (2x3-x1] 

7T 


(2.16c; 


then,  the  quantities  —  can  be  written  in  matrix  form  as 


dxx' 


9x' 


—  0 


2x   /  — 

7T 


2.1 


According  to  equation  2.13  any  contravariant  vector,  say  T.  with  components  T; 
can  be  expressed  in  the  new  (primed)  coordinate  system  as 


N  dx" 


A=l    OX 


or  in  terms  of  components 


T1     =  -v  /— T1  -  0T- 

Tl 


OT 


— T 


(2.18a 


->2' 


OT1 


— T  +  OT" 

7T 


2-l- 


(2.18b) 
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T3'  =   -  x/^-T1  -  OT2  -  2X/^T2       X/^2(T3-T1)  (2.18c) 

4.  Definition  2.4:  Contravariant  Tensor  of  Order  2 

A  set  of  N2  numbers  TA",  where  A  and  ^  =  1.....N  are  said  to  be  the 
components  of  a  second  order  contravariant  tensor  if,  upon  transition  to  another 
coordinate  system,  they  transform  according  to  the  rule 

TA  "    =   E   E  TA"°*     qx  (2.19) 

5.  Definition  2.5:  Covariant  Tensor  of  Order  2 

A  set  of  N2  numbers  UA/1,  where  A  and  //.  =  1,....N  are  said  to  be  the 
components  of  a  second  order  covariant  tensor  if,  upon  transition  to  another 
coordinate  system,  they  transform  according  to  the  rule 

Or,--EE«*!jrr£r  (2-2") 

Similarly,  tensors  of  higher  order  can  also  be  defined.  In  the  general 
case,  it  will  no  longer  be  possible  to  use  different  letters  to  denote  indices.  In  this 
case  indices  with  sub-indices  will  be  utilized,  namely  Ax,  A2.  ....  AN. 

6.  Definition  2.6:  Contravariant  Tensor  of  Order  p 

A  set  of  V  numbers  T  '  p.  where  Xt  =  1,  ...,  N  for  i  =  l,...,p,  are  said  to  be 
the  components  of  a  p-th  order  contravariant  tensor  if.  upon  transition  to 
another  coordinate  system,  they  transform  according  to  the  rule 

Tv    v ,  f.  .      £  T>,    '^  . ,/»*£  (2.21) 

i,    i  .'p- 1  i)x  '  8x  p 

7.  Definition  2.7:  Covariant  Tensor  of  Order  p 

A  set  of  \'p  numbers  I  .        >  .  where  A        l N  for  i       l p.  are  said  to  be 

the  components  of  a  p-th  order  covariant   tensor  if.  upon     transition  to  another 
coordinate  svstem.  thev  transform  according  to  the  rule 


N  ..    p 


dx  ox 


«V      v-    S    ■••     HV     ».fv   ■■■   fv  (2-22) 

A,=  l  A  =i  dx  9x   p 
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8.     Definition  2.8:  Mixed  Tensor  of  Order  p 

A  set  of  Np  numbers  S  l        qA  A    where  A,  =  l N  for  i  =  1 p.  are  said 

to  be  the  components  of  a  p-th  order  mixed  tensor,  with  q  contravariant  and  (p- 
q)  covariant  indices  if,  upon  transition  to  another  coordinate  system,  they 
transform  according  to  the  rule 

=  E  ■••   Ss1     V,-inr-    \     *-     •  •  • -^  (2.23) 

x%        x^i  q+l      p9x '        ax'  ax  <+l        3x  p  v      ; 

We  have  already  seen  two  examples  of  mixed  tensors.     The  first  is  the 
Kronecker  delta 

(1    whenA  -a  ,_  _  „» 

2.24 
0    whenA  #  //  V  ; 

To  see  that  the  Kronecker  delta  is  in  fact  a  mixed  tensor  we  must  test  to  see  if  it 
transforms  according  to  the  rule  given  in  equation  (2.23).  We  must  prove  that 
the  following  relation  is  true 

'''.■=SE^'.  (2-25) 

A=l    M=l      OX        OX 

We  begin  with  the  right  hand  side  of  (2.25) 

(2.26a) 


2.26b 


N       N 

E  E 

dxx' 
9xA 

9x" 

bxM' 

*v 

N 

E 

A-  1 

9xv 
9xA 

N 

9x" 
Bx"' 

= 

N 

E 

A  =  1 

9xA' 
9xA 

3xA 

9x"' 

dx^ 


A' 


-  6*  „-  (2.26d) 

Therefore,  we  conclude  that  relation  (2.25)  holds  and  so  the  Kronecker  delta  is  in 
fact  a  second  order  iensor  of  mixed  character. 
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The  other  mixed  character  tensor  that  we  have  already  worked  with  is 

the  one  appearing  in  the  formulae  of  transformation  (definitions  2.2  through  2.8). 

that  is  the  the  partial  derivative  of  the  new  coordinates  with  respect  to  the  old 

(and  the  old  with  respect  to  the  new).    We  will  not  prove  that  this  is  in  fact  a 

tensor  of  the  type  stated  although  the  proof  is  straight  forward.  It  is  instructive 

3xA  8xA 

to  note,  however,  that  the  two  quantities  — —  and are  inverses  of  each 

8xA  9x^ 

other. 

As  a  final  note,  vectors  are  tensors  of  order  one.  Also  scalars  are 
considered  to  be  tensors  of  order  zero.  They  are  sometimes  called  invariants 
since  their  representation  is  independent  of  the  coordinate  systems  used. 

C.    NOTATIONAL  CONVENTIONS 

There  are  two  widely  accepted  conventions  that  simplify  notation  and 
unquestionably  save  much  writing.  The  first  is  known  as  the  summation 
convention.  Historically,  it  was  first  used  by  Einstein.  He  noticed  that  in  almost 
all  cases  there  is  really  no  need  to  explicitly  write  summation  symbols. 
Summation  can  be  implied  whenever  an  index  is  repeated  in  an  expression,  once 
as  a  superscript  and  once  as  a  subscript.  The  repeated  index  is  allowed  to  take 
on  all  permissible  values  and  the  resulting  terms  are  summed  together.  This  type 
of  index  is  often  referred  to  as  a  dummy  index.  But  what  are  permissible  values 
for  the  index?  This  question  leads  to  the  second  convention.  Normally,  the  range 
of  the  index  will  not  be  explicitly  stated.  By  convention  it  is  understood  that  all 
greek  subscripts  and  superscripts  appearing  in  an  expression  will  take  on  all 
values  from  1  to  N.  where  N  is  the  dimension  of  the  vector  space  in  which  we  are 
working.  In  later  chapters  we  will  find  it  more  convenient  to  allow  indices  to  run 
from  0  To  N.  The  dimensionality  of  the  vector  space  will  thus  be  N-f-1.  An 
additional  convention  which  we  shall  find  useful  is  To  reserve  latin  indices  to 
indicate  ThaT  we  are  dealing  with  a  particular  component  of  a  quantity.  In  most 
books  this  is  indicated  by  surrounding  the  particular  index  with  parenthesis. 
However,  we  will  reserve  parenthesis  to  indicate  exponentiation.  The  conventions 
adopted   here   will   be   used   throughout   the  sequel.     In   exceptional   cases,   where 
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some  deviation  from  them  is  required,  we  will  explicitly  state  the  meaning  of  the 
notation.  As  an  example  of  the  use  of  these  conventions  consider  the  expression 

YA  =    £  HvT*      for  A  =  1,  ...,  N  (2.27) 

It  can  be  written  more  succinctly  as: 

YA  =  HVT*  (2.28) 

Another  convention  which  has  already  been  used  is  now  formally  introduced. 
Every  tensor  quantity  will  be  given  a  distinct  base  letter.  Upon  a  change  of 
coordinates,  the  base  letter  will  be  maintained  in  order  to  indicate  that  the 
quantity  has  not  been  modified,  only  the  representation  has  changed.  The 
coordinate  system  used  is  indicated  by  the  sub-  or  superscript  used  to  index  the 
components  of  the  quantity.  We  therefore,  will  refer  to  different  coordinate 
systems  simply  by  the  index  letter  used  to  indicate  the  components.  For 
example,  a  vector  T  has  components  TA  in  the  (A)  set  of  coordinates,  while  it  has 
components  TA '  in  the  (A')  coordinate  system.  We  note  that  using  this 
representation,  scalars  appear  identical  in  all  coordinate  systems,  which  is 
desirable. 

1.      Example  2.2 

The  following  example  serves  not  only  as  an  illustration  of  the 
conventions  presented  in  this  section,  but  also  as  a  concrete  (and  presumably  a 
somewhat  familiar)  illustration  of  the  two  types  of  vector.  Consider  an  invariant 
function  of  the  coordinates,  f  =  f(x'.  x2.  xs).  The  differential  of  this  function  is  given 
by 

8f   ,  i         Bf    ,  2         Bf 


dx'   ■   --r-dx"  (2.29 

dx'  dx"  dx" 

We  can  consider  dx-*  to  be  the  components  of  a  contravariant  vector  representing 
an  infinitesimal  displacement  expressed  according  to  some  basis  A  =  {a,  a2  *u- 
We  can  write  this  as: 
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dx=    X)  dxAaA  (2.30) 

We   have   already   shown   that    components   of  a   contravariant   vector   transform 
according  to 

dxy  =  -^-dxA  (2.31) 

dxx  K         ' 

The  gradient  is  also  a  vector  whose  components  are  given  by: 

Vf,  =  -~  (2.32) 

Upon  a  change  of  coordinates  the  values  of  the  new  components,  yfA-,  can  be 
deduced  by  application  of  the  chain  rule; 

It   is  clear  that  the  components  of  the  gradient  transform  according  to  the  rule 
given   in  equation    (2.14).     The  gradient   must,  therefore,  be  considered  to  be  a 
covariant  vector. 
2.     Example  2.3 

Although  it  has  already  been  stated  (section  2.1)  that  linear  functionals 
can  be  considered  to  be  covariant  vectors,  this  fact  has  not  been  proven.  In  this 
example  we  will  show  that  any  linear  functional,  say  H.  with  components  HA  that 
transforms  an  arbitrary  contravariant  vector,  say  T.  (according  to  equation 
(2.7b))  to  yield  an  invariant,  satisfies  the  definition  of  a  covariant  vector 
(equation  (2.14)).  Equation  (2.7b).  which  defines  a  vector  inner  product  is 
repeated  here  for  convenience. 

H(T)  =   H,TA  (2.34a) 

The    quantities    T';    are    the    components    of   an    arbitrary    contravariant    vector. 
Because  of  the  assumed  invariance  we  may  write 

HAT  =  HAT'  (2.34b) 

From  the  definition  of  a  contravariant  vector,  (equation  (2.13)) 
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Ty  =  TA-^—  (2.34c) 

Equation  (2.34b)  becomes 

HATA  =  HA-TA  -^—  (2.34d) 

3x 

Rearranging  yields, 

(H>  -  ha-  4nr)  T"  =  °  (2-34e) 

ox 
Equation  (2.34e)  implies 


3x_A' 
8x' 


H,=  Hr^-  (2.34f) 


since  the  contravariant  vector,  T,  was  arbitrary.    A  simple  change  of  variables  in 
this  last  expression  yields 

H,.=  HA-|4  (2-34g) 

ox 

which   is  identical  to  the  equation  defining  a  covariant  vector.    (2.14).     We  are 
thus  justified  in  calling  linear  functionals  covariant  vectors. 

D.    BILINEAR  AND  MULTILINEAR  FORMS 

We  next  consider  higher  order  functionals.  In  particular  we  will  start  with 
the  bilinear  form  or  bilinear  functional. 

1.      Definition  2.9:  Bilinear  Functional 

A  mapping  f  of  a  pair  of  vectors,  say  T  e  UN  and  S  e  VM  into  a  field  of 
scalars.  is  a  bilinear  functional  or  bilinear  form  if  f(T,S)  is  a  linear  function  of 
T  and  S  taken  independently.  Wre  will  only  consider  cases  when  N  =  M  and 
UN   -  VN. 

Once  again  choosing  an  arbitrary  basis  A   =   {a, aN}  we  may  express 

two  arbitrary  contravariant  vectors  as  linear  combinations  of  these  basis  vectors. 

T  =   TAaA,  and  S  -  S"a„ 

The  bilinear  functional  f  can  then  be  written 
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f(S.T)  =  f(S*aM,TAaA)  (2.35a) 

=  S"f(aM,TAaA)  (2.35b) 

=  S*TAf(aM,aA)  (2.35c) 

The  bilinear  form  is  thus  completely  determined  by  the  N2  quantities  f(aM,aA).  We 
will  write  these  components  off  as  fMA.  Using  this  shorthand,  equation  (2.35c)  can 
be  written  as 

f(S,T)  =  f„AS"TA  (2.35d) 

In  matrix  notation  the  bilinear  form  takes  on  the  familiar  appearance 

f(S.T)  =  STFT  (2.36) 

where  F  =  !fMA] 

If  the  two  vectors  S  and  T  are  equal  then  the  bilinear  form  reduces  to  the 
well  known  quadratic  form 

f(T,T)  -  f^T*TA  (2.37a) 

or  in  matrix  notation 

f(T.T)  -   TTFT  (2.37b) 

We     will  be  interested  in  the  behaviour  of    the  components.  f„A.  of  the 
bilinear   form.   f.   upon   transition   from  one   coordinate  system   to   another.     We 
establish  their  tensor  character  in  the  following  example. 
2.      Example  2.4 

In  Example  2.3  we  showed  that  a  linear  functional  satisfied  the  definition 
of  a  covariant  vector.  Here  we  will  show  that  a  bilinear  functional  satisfies  the 
definition  of  a  covariant  tensor  of  second  order.  It  is  necessary  for  the  discussion 
that  follows  in  later  chapters  to  establish  the  tensor  character  of  bilinear 
functionals.  Since  the  bilinear  functional  yields  an  invariant  (scalar),  we  may 
write 

f,A^TA  =  f„-rS"'Tv  (2.38a) 
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=  VA-s^Tir-r  (2-38b) 

dxM  3x 

=  ^—  ^ir-  f„-rS"TA  (2.38c) 

ax"    8xA    M  v        ; 

Rearranging  yields 

(f^77  77Vr)S^=0  .  (2.38d) 

ox'      ox 

Since  the  contravariant  vectors  S  and  T  are  arbitrary  the  quantity  inside  the 
parenthesis  must  vanish.    This  yields  the  relation 

Vr  =^Tirirr  (2.38e) 

oxM      ox 

This  last  equation  is  identical  to  the  definition  of  a  second  order  covariant  tensor 
(equation  (2.20).)  We  have  thus  proven  that  bilinear  functionals  are  covariant 
tensors  of  order  2. 

In  general  we  can  have  m-linear  functionals  which  map  m  contravariant 
vectors.  T(l),  ....  T(m),  into  a  scalar  and  are  linear  functions  of  each  of  the  m  input 
vectors  taken  separately.  They  can  be  written  as 

f(T(l) T(m))=T"I(l)      ••    TAm(m)fv    .  ^  (2.39) 

Using  identical  arguments  to  the  ones  presented  for  linear  and  bilinear 
functionals.  we  can  show  that  multilinear  functionals  are  also  covariant  tensors. 

E.     TENSOR  OPERATIONS 

There  are  a  few  tensor  operations  that  will  be  of  considerable  importance  in 
later  discussions.  Although  some  have  already  been  used  we  will  formally  define 
them  here  before  proceeding.  Only  those  operations  that  will  be  used  in  the 
sequel  are  presented.  Others  are  possible  and  are  discussed  aT  length  in  the 
references  (see  for  example  [Refs.  1,5,6.7].) 
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1.      Definition  2.10:  Tensor  Outer  Product 
Given  the  components  of  two  tensors 


'"p     ]  n 


'1  ^m 


A       ,  A 

and  S  n  +  1         p, 


the  N(p+q)  numbers  RAl      x\i .. .  M   given  by 


"l  "q    ~  "l  MmD  Mm+1         "<, 


:2.40] 


(2.41) 


are  components  of  a  tensor  of  order  p  +  q.    The  operation  implied  in  the  above  is 
known  as  the  tensor  outer  product  or  simply  the  tensor  product. 

2.     Example  2.5 

As    an   example    of  the   tensor   product    operation    consider   two   vectors 
defined  as: 


T  =   [TA:        and        S  =   [S*] 
The  tensor  product  is  given  by  (equation  (2.41)) 


(2.42) 


(2.43) 


In  this  case  the  components  of  the  tensor  product  can  be  arranged  as  a  matrix  of 
N2  numbers.  For  simplicity  consider  the  case  when  N=3.    In  matrix  notation 


or 


R  =   TS1 


RA„    =    TAS"> 


;2.44a) 


[2.44b) 


T'lci  'T'Q'  T'S* 
T-S1  T2S2  T2S; 
T'S1    T'S-    T3S' 


Definition  2.11:  Contraction 


!.44c 


The  operation  of  setting  two  indices,  one  lower  and  the  other  upper, 
equal  and  summing  the  result  is  known  as  contraction.  The  result  is  a  tensor 
which  has  the  character  indicated  by  the  remaining  indices.  The  contraction  of  a 
tensor  of  order  p  +  1  over  two  indices  results  in  a  tensor  of  order    p  -  1. 
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4.     Example  2.6 

We  may  contract  the  tensor 


H^ 


H  j  H  2  h  j 

H2  u2  u2 

i  ±1  2  rl  3 

H3  ttS  tt3 

1  «   2  "   3 


(2.45) 


over  the  indices  A  and  /*,  resulting  in 
H\  =  H\  +  H22  +  Hss 


(2.46a) 


=  Trace|H; 


(2.46b) 

a  scalar  which  represents  the  trace  of  the  matrix. 

We  can  consider  a  higher  order  example.  Assuming  a  three  dimensional 
vector  space,  consider  contracting  a  tensor  UA^7  over  the  indices  A  and  p.  The 
result  will  be  a  vector  whose  components,  IP,  are  given  by 

U1  =  UV  +  UV  +  U3,1  (2.47a) 


U2  =  U1 


u2,2  +  us  2 


Us  =  UV  +  UV  -  u3 


5.     Definition  2.12:  Tensor  Inner  Product 


;2.47b) 


(2.47c) 


Suppose  that  after  we  form  the  outer  product  of  two  arbitrary  tensors  T 
and  S.  with  components  T  '        '        .  u     and  S  nTl        pu  „  .  we  set  the  indices  A. 

and    |ij    equal.     This    implies     a    contraction    operation.     The    outer    product    is 
written  (according  to  equation  (2.41))  as 


RA' 


T*. 


sA- 


'2.48a) 


The  contraction  operation  yields.  (  definition  2.11) 


RA' 


A:     ,A..l 


j     l/j-J 


1-   ,'   j-1 


2.4S1> 


It  can  be  shown  that  the  object.  R.  given  in  the  last  equation  is  a  tensor  of  the 
character  shown.  If  the  original  tensors.  T  and  S,  were  of  order  m-f-n  and  p  +  q- 
(m+n)   respectively,  then  the  resulting  tensor.  R.  will  be  of  order  (p  +  q-2).   The 
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total  operation,  consisting  of  contracting  the  result  of  a  tensor  product  over  one 
(or  more)  pairs  of  indices,  each  member  of  the  pair  having  come  from  a  different 
tensor,  is  known  as  a  tensor  inner  product,  or  simply  inner  product.  This 
operation  is  also  sometimes  referred  to  as  transvection  (see  for  example  [Ref. 
5].) 

This   operation   has    been    used   previously    in    the    discussion   of  linear, 
bilinear   and   multilinear   functionals,    (see    equations    (2.7b),    (2.35d)    and    (2.39) 
respectively.)    It  will  be  particularly  valuable  in  future  discussions. 
6.      Example  2.7 

Some  insight  into  the  inner  product  operation  may  be  gained  by 
explicitly  performing  the  two  steps  described  above  for  the  simple  case  of  the 
linear  functional  (equation  (2.7b).) 


y  =  H,T; 


f2.49l 


We  can  first  form  the  outer  product  by  replacing  one  of  the  A  indices  appearing 
on  the  right  hand  side  with  a  different  letter,  say  fi .    We  may  then  write: 


!H,T> 


H,T'  H,T2  H,T3 
H2T'  H2T2  H2T3 
H,T]    H,T2    H,T3 


'2.501 


The  result   is  now  contracted  by  equating  A   and  //    and  performing  the  implied 
summation. 


y  "  H,TA 


2.51a 


H,T]  -  H-.T- 


H.r- 


2.51b 


It  is  often  useful  to  perform  both  these  steps  (tensor  outer  product  and 
contraction)  when  calculating  an  inner  product,  particularly  when  dealing  with 
higher  order  tensors.  It  may  otherwise  be  difficult  to  keep  track  of  how  terms  are 
combined  or  even  which  terms  will  be  present  in  the  final  expression. 
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III.    NONLINEAR  SYSTEMS  THEORY 

In  this  chapter  several  models  of  nonlinear  systems  are  described.  The 
discussion  begins  with  an  overview  of  the  classical  continuous  Volterra  series. 
Several  interesting  aspects  of  the  series  are  examined.  Next,  a  discrete-time 
version  of  the  Volterra  series  is  introduced  and  is  used  to  develop  an  equivalent 
tensor  form.  Finally,  a  new  discrete-time,  nonlinear  tensor  model  is  presented 
and  its  relationship  to  the  Volterra  series  is  discussed. 

A.    CONTINUOUS  NONLINEAR  SYSTEMS  MODELS 

Traditionally,  non-linear  systems  have  been  modelled  using  the  continuous 
Volterra  series  expansion  [Ref.  26:  p.  1559].  In  its  most  general  form  this  can  be 
written  as: 

y(t)  =  h0(t) 

OG 

-  'fb.l{t,Tl)x[rl)dr1 


J   J  h2(t,7-1.r2)x(7-1)x(r2)dr]d7-2 


-    •  •  •  (3.1) 

where  x(t)  is  the  system  input  and  y(t)  is  the  system  output.  The  parameters 
h0(t).  h,(t.r,).  h2(t. r,.r2).  are  known  as  the  Volterra  Kernels.    As  we  will  only 

be  concerned  with  time-invariant  systems,  the  kernels  will  only  be  functions  of 
the  time  difference,  (t-r).  and  not  the  actual  time.  The  kernels  then  become: 
h0,  hift-rj).  h2(t-r1,t-T2).  ■■-.  Simple  changes  of  variables  then  allow  the  series. 
(3.1),  to  be  written  in  the  form 
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y(t)  =  h0 

+  J  hi(riWt-ri)dr, 

—  00 

00      00 

+  J  J  h2(r  !,t 2)x(t-T  Jxit-r  2)dr  tdT  2 

—  00—00 

+    •••  (3.2) 

We  note  several  things  about  the  expansion.  First,  an  infinite  number  of 
terms  are  required  to  represent  the  most  general  case.  Second,  the  kernels 
h0.  h,(7-];  h2(r1,r2),  ■  '  correspond  to  the  constant,  linear,  quadratic,  ...  terms  of 
the  expansion,  respectively.  The  familiar  linear  system  appears  as  a  special  case 
of  this  more  general  expansion  (ie;  the  case  when  all  kernels  except  h^r,)  are 
zero.)    The  third  term: 

OC      00 

J   J  h2(7",,r2)x(t-7-1)x(t-r2)dr1dr;,  (3.3) 

—  oo—  oc 

is  a  bilinear  term,  that  is  it  is  linear  in  each  variable  x(t-rj)  and  x(t-r2)  taken 
independently  (ie;  assume  the  other  variable  is  a  constant.)  In  general,  the 
(i+l)-st  term  is  i-linear.  It  is  linear  in  each  of  the  i  variables 
x(t— rj).  x(t— t2),  •••,  x(t-ri)  taken  one  at  a  time.  Lastly,  notice  that  the  Yolterra 
expansion  is  not  orthogonal  in  the  sense  that  the  identification  of  the  n-th  order 
kernel  depends  on  the  values  of  all  the  other  kernels.  They  cannot  be  identified 
independently. 

The  non-linear  model  represented  by  equation  (3.1)  can  be  visualized  as 
shown  in  Figure  3.1.  As  can  be  seen  it  corresponds  to  a  parallel  connection  of 
subsystems  of  increasing  non-linearity.  Each  of  the  subsystems  is  homogeneous 
(except  for  the  zeroth  order  subsystem)  in  the  sense  that  increasing  (multiplying) 
the  input  by  a  factor  k  results  in  the  output  of  the  p-th  subsystem  being 
increased  by  a  factor  kp.  This  can  easily  be  understood  by  examining  the 
expression  for  the  p-th  term. 
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OC  00 


J  •■  J\(t- !,...,  rp)x(t-r,)       ••    x(t-rp)drj    ■■•    drp  (3.4) 

—  00  —00 

Replacing  x(t)  by  kx(t)  yields 

00  00 

j  ■  ■  ■  Jhp(r„  ...,  rp)kx(t-ri)    ■•    kx(t-rp)dr,    ■■■    drp  (3.5a) 


-00  —00 


=   k"J  •  •  •  Jhp(rl5  ...,  rp)x(t-r,)    ■  •  ■    x(t-rp)dr,    •  ■  ■    drp  (3.5b) 

—  00  —00 

The  presence  of  a  constant  term  in  the  Volterra  expansion  should  not  be 
unexpected.    Consider,  for  example,  a  system  whose  response  is  given  by: 

y{t)=x(t)+h  (3.6) 

It  is  easily  shown  that  this  system  does  not  obey  the  principle  of  superposition 
and  so  cannot  be  considered  linear.  This  necessitates  the  inclusion  of  a  constant 
term  in  the  Volterra  expansion  in  order  to  handle  the  general  case.  The  usual 
procedure  adopted  in  linear  analysis,  if  a  constant  term  appears,  is  to  define  a 
new  output  function  which  is  the  actual  output  less  the  constant  term.  This  new 
output  function  is  then  identified  in  the  usual  fashion.  The  constant  term  must 
be  separately  identified.  If  we  admit  that  the  system  is  non-linear  then  this  will 
no  longer  be  necessary. 

There  are  many  other  aspects  of  continuous-time  nonlinear  system  modelling 
that  have  not  been  discussed.  In  the  next  section  discrete-time  nonlinear  systems 
are  introduced.  Many  of  the  comments  that  will  be  made  there  are  equally 
applicable  to  continuous-time  nonlinear  systems.  However,  we  make  them  in  the 
context  of  discrete-time,  since  that  will  make  them  more  applicable  to  the  sequel. 

B.    DISCRETE-TIME  NONLINEAR  SYSTEM  MODELS 

A  discrete  version  of  the  time-invariant  Volterra  expansion  can  be  written  as: 


y(k) 


EMA^k-A,) 
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Figure  3.1:  Nonlinear  Yolterra  Series  Model 
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-    E     E  hatA^AJxtk  -A,)x(k-A2) 


(3.7) 

where  we  have  assumed  that  the  system  is  causal. 

For  a  large  class  of  systems  the  summations  can  be  truncated  to  N-t-1  terms. 
Certainly  truncation  is  required  for  computer  implementation.  This  truncation 
implies  that  only  finite  memory  nonlinear  systems  can  be  represented.  Equation 
(3.7)  is  an  extension  of  the  linear  Moving  Average  (MA)  type  model.  In  fact  the 
linear  model  is  a  special  case  of  equation  (3.7).  This  expansion  is  non-recursive  in 
nature,  that  is.  it  expands  the  present  output  only  in  terms  of  the  present  and 
past  input.    Past  outputs  are  not  used. 

The  representation  of  the  kernels  in  both  the  continuous  and  discrete  forms 
of  the  Volterra  expansion  is  not  unique.  There  are.  however,  several  special  forms 
which  are  important.  Consider  a  second  order  kernel  for  which 
MAi,A2)  =  h2(A2,Ai).  This  kernel  is  symmetric  with  respect  to  the  two  parameters  A, 
and  A2.  It  turns  out  that  the  kernel  can  always  be  symmetrized  with  no  loss  of 
generality.  For  the  p-th  order  kernel  the  procedure  for  obtaining  the  symmetric 
kernel  from  an  asymmetric  one  is  given  by  [Ref.  22]: 

tsym(Aj A2)  =  — r£hp(A*(0'  ■■■■.  A„(p))  (3.8) 

n-   r() 

where  the  summation  is  over  all  n!  possible  permutations  of  the  p  A  's.  Although 
the  symmetric  form  may  contain  more  terms  than  an  asymmetric  form  it  is  of 
importance  because  it  is  unique  [Ref.  9:  p.  43].  There  can  be  many  equivalent 
asymmetric  forms  of  the  kernel  which  all  lead  to  the  same  symmetric  kernel 
(through  equation  (3.8)).  The  symmetric  kernel  thus  provides  a  standard  form 
which  can  be  used  as  a  reference. 

Other  forms  of  interest  are  also  possible.  The  symmetry  of  the  symmetric 
kernel  implies  redundancy.  This  redundancy  can  be  eliminated  by  use  of  a 
triangular  kernel.     Consider  a  kernel  defined  so  that  h,n(A, Ar)       0  whenever 
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A9>A,  for  s  <  t.  The  domain  of  a  second  order  triangular  kernel  is  illustrated  in 
Figure  3.2a.  For  comparison,  the  equivalent  symmetric  kernel's  domain  is 
illustrated  in  Figure  3.2b. 

Using  the  triangular  kernel  the  output  of  the  p-th  order  nonlinear  subsystem 
is  given  by 


N 


>'(k)=    E    E     •■•     Ehtri(A„...,Ap)x(k-A1)      •x(k-Ap)  (3.10) 


v° 


Notice  that  the  limits  of  the  summations  reflect  the  triangular  domain.  This 
implies  that  fewer  terms  are  included  in  the  summations  resulting  in 
computational  savings.  The  triangular  kernel  defined  above,  and  used  in 
equation  (3.10).  is  not  unique.  Other  triangular  kernels  can  be  formed  by 
choosing  alternate  triangular  domains.  For  example,  the  domain  illustrated  in 
Figure  3.3  could  equivalently  have  been  used.  This  choice  corresponds  to  setting 
h  tn(Ai,   .-,  Ap)  =  0  whenever  A5  >  At,  for  s  >  t. 

The  output  of  nonrecursive  models  of  the  type  presented  in  equation  (3.7)  is 
stable  if  the  input  is  bounded  and  if  the  series  is  truncated  to  a  finite  number  of 
summations.  Consider  an  input  x(n)  which  is  bounded  to  be  less  than  some 
constant  M.  If  the  series  is  truncated  to  p-fl  terms  then,  in  the  worst  case  the 
output  will  be 

y(nKh+M£      MA,)      -  M2  £     £  i  h,(A,.A2)  ;    + 

A,=  0  A,=  0   A2=0 

•      -  MPE    •••     E     M*. Ap)  (3.11) 

So.  as  long  as  the  summations  are  bounded  (which  will  generally  be  the  case),  the 
output  will  always  remain  finite.  This  guaranteed  stability  makes  MA  type 
models  very  attractive.  As  mentioned  earlier,  their  shortcoming  is  their  inability 
to  accurately  model  infinite  memory  systems  without  using  a  large  number  of 
terms. 
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a.    Second  Order  Triangular  Kernel  Domain 


b.    Second  Order  Symmetric  Kernel  Domain 


Figure  3.2:  Triangular  and  Symmetric  Yolterra  Kernel  Domains 
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1.     Tensor  Formulation 

In  order  to  adopt  a  tensor  formulation  of  the  problem  we  notice  that 
equation  (3.7)  can  be  considered  to  consist  of  a  series  of  increasing  order 
functionals.  As  has  been  shown,  these  functionals  can  be  expressed  as  tensors. 
We  can  therefore  rewrite  equation  (3.7)  as  a  tensor  equation. 

y(k)  =  H 


+  HAixAl 


+  Hjj  A  x  'x 


+   •  •  •  (3.12) 

where  we  have  defined  the  contravariant  input  vector  as 

x  =  ix(k)  x(k-l)    ■••    x(k-A)    ■••    x(k-N)]T  (3.13a) 


=   ;xA,T 


(3.13b) 


This   choice   for   x   has   the   effect    of  truncating  the   series  to   N+l   terms.     The 
symbols  H,  HA  .  HA  A  .  represent  the  components  of  covariant  tensors  of  order 

0.  1,  2,  etc..  respectively. 

Examining  equation  (3.12)  we  make  particular  note  of  the  following  two 
aspects: 

(1)  The  dimension  of  the  vector  space  is  related  to  the  memory  order  of  the 
system.  The  vector  x  has  N  +  l  components  implying  that  the  nonlinear 
system  contains  no  more  than  N  delays.  This  fact  is  explicit  in  the  way  the 
input  vectors  have  been  defined. 

(2)  The  nonlinearity  of  the  expansion  is  provided  implicitly  by  the  tensor  outer 
product  operation.  It  is  the  outer  product  of  the  vector  x  with  itself  that 
makes  the  higher  order  (2  and  up)  terms  nonlinear. 
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Figure  3.3:  Alternate  Second  Order  Triangular  Yolterra  Kernel  Domain 
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These    observations    provoke    speculation    about    the    possibility    of    an 
alternate  formulation  where  the  roles  played  by  the  dimensionality  of  the  vector 
space  and  the  outer  product  operation  are  interchanged. 
2.     Alternate  Tensor  Formulation 

Consider  a  parametric  description  of  a  curve  in  Vp+1,  a  p+1  dimensional 
vector  space,  given  by; 

x°(k)  =   [x(k))<°> 
x*(k)  =  [x(k)p) 


x(k)  = 


xp(k)  =  [x(k)](p) 


(3-14) 


where  the  superscript  in  parenthesis  indicates  exponentiation.  Given  any  x(k). 
the  components  x"*(k),  A  =  0,...,p  define  a  point  in  Vp+1.  As  x(k)  varies  with  k.  the 
vector  x(k)  will  describe  a  curve  in  Vp+1.  We  define  the  components  of  another 
vector  in  Vp+1  as 


xAi(k-i)  =  ;x(k-i)](Ai) 


Similarlv.  the  N-th  such  vector  can  be  defined  as. 


x"N(k-N)  =   |x(k-N)!(AN) 


(3.15) 


(3.16) 


The  vectors  in  equations  (3.14).  (3.15).  and  (3.16)  will  be  referred  to  as 
observation  vectors.  Although,  at  this  point  they  only  depend  on  past  and 
present  system  input  values,  later,  in  Section  D,  they  will  be  generalized  to 
include  past  outputs  as  well.  The  input  and  output  measurements  represent  the 
system  observations,  hence  the  name. 

Using  the  theory  developed  in  Chapter  2.  concerning  multilinear 
functionals,  the  following  mathematical  relation  is  proposed  as  a  model  of  a  finite 
order,  finite  memorv  nonlinear  svstem: 


y(k)  -  x  "(k) 


x  N(k-N)HAr 


(3.17) 
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where  HA  .  A  is  an  (p-l-l)-order  covariant  tensor  representing  a  (p+l)-linear 
functional.  This  covariant  tensor  performs  a  role  similar  to  that  of  a  Volterra 
kernel. 

This  model  is  equivalent  to  the  p+1  term  Volterra  series  (equation  (3.7) 
or  (3.12)),  although  it  does  contain  many  additional  terms.  It  has  the  advantage 
of  notational  simplicity.  It  replaces  a  p+1  term  series  with  a  single  term  and 
requires  the  specification  of  only  a  single  composite  kernel,    HA  . .   A  .    As  will  be 

shown  in  Section  D,  equation  (3.17)  is  considerably  more  general  than  equation 
(3.12).  It  will  be  shown  that  Wiener  type  models  can  be  obtained  from  equation 
(3.17)  by  a  simple  coordinate  transformation.  Other  choices  of  observation 
vectors  yield  Autoregressive  (AR)  or  even  Autoregressive-Moving  Average 
(ARMA)  type  models. 

In  order  to  illustrate  the  correspondence  of  equation  (3.17)  to  the 
standard  Volterra  type  series  expansion  of  equation  (3.12)  we  present  a  simple 
example. 

3.     Example  3.1 

Consider  the  truncated  Volterra  series  expansion  corresponding  to 
equation  (3.12); 


y(k)  =   H  +  HA/ 1  -  HAiVV2 
where  Aj  =  0.1    and  A2  =  0,1.  and 


(3.18) 


x(k) 

x(k-r 


We  may  explicitly  write  out  all  the  terms  implied  in  equation  (3.18).    This  yields 
y(k)  =   H 

-  H0x(k)  +  H,x(k     1) 


+  H00x(k)x(k)  ~  H0,x(k)x(k-i; 


+  H,0x(k-l)x(k)  +  H,,x(k-l)x(k-l) 


(3.19) 
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Now  consider  a  model  of  similar  order  given  by  equation  (3.17). 

y(k)=  Hv/°(k)xAl(k-l)  (3.20) 

where  A0,  Aj  =  0,1,2    and 

xA°(k)  =  [x(k)](Ao) 

xAl(k-l)=   |x(k-l)](Al) 

The    tensor   product,    x  °(k)x  '(k-1),     appearing    in   equation    (3.20)    results   in    a 
contravariant  tensor  of  second  order.  Its  elements  are 


|xAo(k)xV-l)]  = 


1  x(k-l)  x(2»(k-l) 

x(k)       x(k)x(k-l)       x(k)x(2,(k-l) 

x'2)(k)    x(2>(k)x(k-l)    x(2)(k)x'2'(k-l) 


(3.21) 


We  notice  that  all  the  elements  on  or  above  the  southwest-  northeast  diagonal 
are  included  in  equation  (3.19).  the  terms  below  this  diagonal  are  not.  This  new 
form  has-  not  discarded  any  terms  present  in  the  latter  version  and  so  cannot 
represent  a  loss  of  generality.  The  extra  terms  that  are  included  in  this  new  form 
do  not  pose  any  significant  problem.  If  the  system  does  not  contain  terms 
involving  these  particular  elements  then  the  corresponding  term  in  the  kernel  will 
go  to  zero  during  the  identification  process.  Certainly  in  some  classes  of  systems 
these  additional  elements  will  be  important. 

C.    DISCRETE  NONLINEAR  SYSTEM  IDENTIFICATION 

In  Section  B  a  tensor  formulation  of  the  discrete  nonlinear  modelling  problem 
was  introduced.  In  this  section  methods  for  kernel  identification  are  presented. 
All  the  methods  that  are  discussed  are  statistical  in  nature  and  utilize  a  least- 
squares  approach  of  parameter  estimation.  It  will  be  shown  that  familiar 
methods  used  in  linear  systems  modelling  can  be  extended  to  handle  the 
nonlinear  case.  In  the  first  section  a  tensor  equivalent  of  the  normal  equations. 
which  can  be  solved  for  the  unknown  system  parameters,  is  derived.  Several 
recursive  (in  time)  solutions  to  the  problem  are  then  presented.  Finally, 
simulation  results  are  offered  to  illustrate  the  effectiveness  of  the  algorithms. 
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1.      Derivation  of  Normal  Equations 

In  Section  B  (equation   (3.17)).  the  following  tensor  model  for  nonlinear 
systems  was  proposed: 


y(k)  =  x°(k)    ■••   xN(k-N)Hv     An  (3.22) 

where  y(k)  represents  an  estimate  of  the  system  output  and  the  x  '(k— i)  's  are 
components  of  contravariant  observation  vectors,  defined  in  equations  (3.14) 
through  (3.16). 

Consider  the  analysis  model  of  Figure  3.4.  This  diagram  represents  a 
conceptualization  of  the  system  identification  process.  The  assumption  is  made 
that  the  unknown  system  can  be  represented  by  an  equation  identical  to  the 
model  equation,  (3.22),  where  the  system  parameters  HA  ...x    are  unknown.  The 

parameters  of  the  model  are  adjusted  to  best  match  the  actual  system 
parameters.  A  convenient  measure  of  how  well  the  model  represents  the  actual 
system  is  the  mean-square  error  or  MSE.  The  MSE  is  a  quadratic  function  of  the 
model  parameters  which  implies  that  there  exists  a  unique  minimum,  or  optimal 
solution,  to  the  problem.  Minimizing  the  MSE  yields  a  linear  set  of  simultaneous 
equations  which  can  be  solved  for  the  unknown  model  parameters.  In  addition, 
the  quadratic  nature  of  the  MSE  allows  steepest  descent  type,  adaptive 
algorithms  to  be  used. 

The  error  signal,  e(k).  defined  as  the  difference  between  the  actual 
nonlinear  system  output.  y(k).    and  the  output  of  the  model.  y(k),  is  given  by 

e(k)=y(k)-y(k)  (3.23) 

The  mean-square  value  of  this  error  signal  is  given  by 

EJe-(k)     =    EJy(k)  ->(k)  -  (3.24a) 


|         E/ryjKj-jiKj     . 


A A, 


-N)HAq      xA  (3.24b) 
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Figure  3.4:  System  Analysis  Model 
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where  E{  }  indicates  expectation.    The  optimal  set  of  parameters  can  be  found  by 
minimizing   the     MSE,   E{e2(k)}.   with   respect    to   the   parameters    HA        A  ..     The 

gradient  of  the  MSE  is  formed  by  differentiating  with  respect   to  the  unknown 
parameters. 

9E|e2(k)l 


9HA 


An  '   A* 


=  2E|[y(k)  -  x  «(k)         x  N(k-N)HAo      AN][-xM°(k)  ■  ■  ■  x"N(k-N)]  j  (3.25a) 

Setting  this  equal  to  zero  yields 
2E/;y(k)x'i0(k)...xMN(k-N)  -  xA°(k)...xAN(k-N)x"°(k)...xMN(k-N)HAo       A 

=  0  (3.25b) 

or. 

E|xA°(k)...xAN(k-N)x'lo(k)...xMN(k-N)|HAo.     An  =  E|y(k)xMo(k)...XMN(k-N)j 

(3.25c) 

The  expression.  (3.25c).  is  a  tensor  equivalent  of  the  Wiener-Hopf  or 
normal  equations.  As  will  be  shown  in  Section  D.  these  equations  ran  also  be 
used  to  represent  the  Yule- Walker  equations  if  a  different  choice  is  made  for  the 
observation  vectors,  x(k-i).  The  Term  on  the  left-hand-side  of  equation  (3.25c)  is 
a  nonlinear  extension  of  the  autocorrelation  matrix.  This  contravariant  tensor  of 
order  2(N-fl)  includes  various  higher  order  correlations  as  well  as  the  second 
order  ones  which  arise  in  the  linear  case.  The  expectation  on  the  right-hand-side 
represents  a  cross  correlation  between  the  output  and  various  linear  and 
nonlinear  functions  of  the  input. 

The  minimum  MSE  can  be  determined  by  substituting  (3.25c)  into 
(3.24b).    This  yields 
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E|e2(k)|      =  EJy2(k)|  -  2EJy(k)XA°(k)  •  ■  ■  x^(k-  N)HAq       A  \ 


^HAo      AnH,o      MNE(xA°(k)...x^(k-N)x"0(k)...x'iN(k-N)  (3.26< 


=  E  y2(k)     -  dy(k)xA°(k)  •  •  ■  x'N(k-  N)   HAn       ,  (3.26b) 


Equation  (3.25c)  represents  (p^l)(N  +  1)  equations  in  as  many  unknowns, 
and    so    can    theoretically    be    solved    for    the    unknown    parameters    HA  . . .  A  . 

However,  in  practice  the  number  of  computations  required  to  perform  the  matrix 
inversion  becomes  unwieldly.  An  nxn  matrix  inversion  takes  on  the  order  of  0(n3) 
operations  [Ref.  27:  p.  58],  therefore,  0([p+l]s^N+1^)  operations  will  be  required  in 
the  solution  of  (3.25c).  In  order  to  make  the  task  manageable,  alternate 
algorithms  that  avoid  direct  matrix  inversion,  must  be  employed  to  solve  the 
normal  equations. 

2.      The  Least  Mean  Square  (LMS)  Algorithm 

The  Least  Mean  Square  (LMS)  algorithm  has  successfully  been  employed 
in  the  solution  of  a  variety  of  linear  problems  [Ref.  28].  It  is  a  recursive 
algorithm  based  on  a  gradient,  steepest  descent  type  of  strategy.  The  mean- 
square  error  is  a  hyperparaboloid  which  is  concave  upward.  Steepest  descent 
algorithms  strive  to  descend  down  this  quadratic  surface,  towards  the  minimum, 
by  making  adjustments  to  the  unknown  parameters  proportional  to  the  gradient. 
This  can  be  expressed  mathematically  as 

Hv     As.(k-])  =  HAo      AN(k)-/,vv     AN(k)  (3.27) 

where  Hs  \  (k)  is  the  value  of  the  model  parameters  at  the  k-th  time  instant 
and  fj  is  a  parameter  which  controls  the  convergence  of  the  algorithm.  The 
symbol.    ;A        ^  ,(k).  is  used  to  represent  the  gradient. 

The  actual  value  of  the  gradient  is  given,  according  to  equation  (3.25a). 
by 
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^o       AN(k)=        2E|e(k).xA°(k)  •     V»(k-N)j  (3.28) 

The  LMS   algorithm  approximates  the   MSE  by   the  instantaneous  value  of  the 
error  squared.    The  approximate  value  of  the  gradient  is 

VA0.AN(k)  =   -  2e(k)[xA°(k)  •  •  •  x^k-N)]  (3.29) 

where  the  circumflex  indicates  an  estimated  value.    Using  this  approximation  in 
equation  (3.25)  yields 

H,o      AN(k-l)  =  Hv     ,N(k)  +  2Me(k)xA°(k)  •     -x^(k-N)  (3.30) 

Equation  (3.30)  gives  a  straight  forward  method  of  determining  the 
system  parameters.  It  involves  no  matrix  inversion  nor  does  it  require  the 
correlation  tensor  to  be  known.  These  two  properties  make  the  calculations 
required  at  each  iteration  very  simple,  so  that  it  is  possible  to  perform  them  in 
real  time.  It  can  be  shown  that  the  LMS  algorithm  converges  to  the  optimal 
solution  [Ref.  29:  p.  578]. 

For  the  linear  case,  the  algorithm  will  converge  as  long  as  the  parameter. 
n .  is  chosen  to  satisfy.  [Ref.  29:  p.  578]. 

Cx^-.-J—  (331) 

■A. 


'max 


where  Ama„  is  the  largest  eigenvalue  of  the  correlation  matrix.  Since  the 
correlation  matrix  is  positive  definite,  all  its  eigenvalues  are  greater  than  zero. 
Therefore,  the  trace  of  the  correlation  matrix  will  always  lie  greater  than  the 
largest  eigenvalue1.    The  following  condition  will,  therefore,  ensure  stability. 


0<, 


Tr;  correlation  matrix 
For  the  nonlinear  case  this  translates  to 

l 


(K/i 


£    ■••    EE{xA°(k)...xAN(k    N)xA"(k)...x^(k-.\)  } 
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In    practice    a    value    considerably    smaller    than    the    maximum    permitted    by 
equation  (3.33)  is  used  to  give  slower,  more  accurate  convergence. 
3.     Recursive  Least  Squares  (RLS)  Algorithm 

The  recursive  least  squares,  or  RLS  algorithm,  is  similar  to  the  well 
known  Kalman  filter  except  that  it  is  used  to  estimate  model  parameters  rather 
than  the  system  state.  The  tensor  form  of  the  normal  equations  is  not  suitable 
for  RLS  implementation.  The  elements  of  the  correlation  tensor  must  first  be 
rearranged  in  a  two  dimensional  matrix.  This  can  be  accomplished  by  replacing 
the  tensor  outer  products  appearing  in  equation  (3.25c)  with  matrix  Kronecker 
products.  The  Kronecker  product  of  an  n>m  matrix,  A=  [a>  A  ],  and  a  p*q 
matrix.  B  =  |bM  „  j,  is  an  np  mq  matrix  given  by  [Ref.  30.31], 


^ 


anB    a12B 
a2]B    a22B 


anIB    a^B 


almB 

a2mB 


a_B 


(3.34) 


where  the  symbol    :L    is  used  to  denote  the  Kronecker  product  operation. 

In  order  to  rewrite  the  normal  equations.  (3.25c).  in  this  matrix  form,  the 
covariant  tensor  of  system  parameters,  HA  x  ■  must  be  put  into  a  vector  form 
by  a  lexicographic  reordering  of  the  elements.  The  resulting  parameter  vector  is 
denoted  H.  The  normal  equations  become 


El  x(k)  :<■•■■    Qx(k-  N)   £  x(k) 


&x(k     1)TH 


-   E<v(k)  x(k) 


x(k-N 


(3.35; 


If  an  assumption  of  ergodicin    i.»  mad<   then  the  statistical  averages  in  (3.35)  can 
be  replaced  with  timi    averages. 


lini   

K-DC       K 


£  x(k)xT(k) 

k^  0 


H 


lini  — 

K  -ex     K 


E>(k)X(k) 

k  =  0 


(3.36) 
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where  X(k)  =    x(k)   ®  •  •  ■   £■  x(k    N)  .   For  computational  purposes  equation   (3.36) 
is  approximated  as 

S  X(k)XT(k)  =    £  y(k)X(k)  (3.37) 

k=0  k=0 

As  a  further  notational  convenience,  the  following  two  definitions  are  made, 


XT(0) 
XT(1) 


XK  = 


3f(K] 


(3.38) 


and. 


y(o) 
y(i) 


y(K) 


(3.39) 


Using    these    definitions    the    normal    equations    can    be    written    in    the 
compact  form 

XKTXKH  =  XKTY  (3.40) 

This  last  from  of  the  normal  equations  is  precisely  the  same  as  the  one  used  by 
Goodwin  and  Payne  [Ref.  32:  p.  176]  for  the  derivation  of  the  RLS  algorithm. 
The  derivation  involves  the  use  of  the  matrix  inversion  lemma  |Ref.  33:  p. 
247]  which  replaces  a  matrix  inversion  at  each  iterative  step  by  a  simple  scalar 
division.  It  is  this  simplification  which  yields  the  efficiency  of  the  RLS  algorithm. 
The  RLS  equations  arc-  [Ref.  32:  pp.  176-177]. 


HK^,       HK   -   QK.,(yiK)       X'(K-  i)H 


.41a 


!h     1 


PKX(K-  li 
XT(K-  1)PKX(K-  1 


(3.41b) 
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Pk  +  i=  [I-Qk^Xt(K-1)Pk  (3.41c) 

=  !XI?+1XK+1]-1  (3.41d) 

Equation  (3.41a)  deserves  some  comment.  The  term  in  the  parentheses 
is  normally  known  as  the  innovation.  It  is  a  scalar  which  represents  the  new 
information  gained  in  the  latest  measurement.  The  term  XT(K+1)HK  is  an 
estimate  of  the  current  output  given  the  new  input  measurement  and  the  old 
estimate  of  the  system  parameters.  The  innovation  is  thus  an  error  signal.  The 
term  in  front  of  the  bracket,  QK+i,  is  a  gain  vector.  It  gives  an  indication  of  how 
much  faith  is  being  put  in  the  new  information.  If  the  gain  is  small,  then  the 
new  estimate  will  be  essentially  the  old  estimate.  Conversely,  if  the  gain  is  large 
then  the  new  estimate  will  depend  to  a  great  degree  on  the  new  information. 
When  the  algorithm  has  converged,  the  gain  will  be  close  to  zero.  If  the  system 
parameters  change  for  any  reason,  the  algorithm  will  not  detect  the  change  since 
it  is  ignoring  the  new  ■  information.  In  order  to  circumvent  this  problem  a 
weighting  can  be  applied  to  the  data,  so  that  new  data  is  artificially  emphasized. 
Exponential  and  rectangular  windows  have  been  successfully  employed  for  this 
purpose.  The  time  constant  used  in  the  case  of  the  exponential  window  is  often 
called  a  forgetting  factor  since  it  ensures  that  the  algorithm's  memory  is  finite. 
Use  of  windows  will  not  be  discussed  further  in  this  dissertation.  The  interested 
reader  should  consult  reference  [Ref.  32:  pp.  179-185]. 
4.      Simulation  Results 

In  order  To  investigate  the  validity  of  the  Theoretical  results,  several 
computer  simulations  were  performed.  Two  examples  are  presented,  representing 
Two  different  classes  of  systems.  Many  other  systems  were  also  tested,  however. 
tIk  results  presented  are  Typical  of  those  obtained  from  all  the  simulations.  The 
FORTRAN  programs  written  allowed  nonlinearities  of  up  To  fourth  order,  but 
were  limited  To  systems  involving  only  zero  or  one  delay. 

The  first  example  was  chosen  to  correspond  exactly  to  the  model 
equation    (3.17).     The  system  was  excited  by  white,  zero  mean,     uniform  noise. 
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The  problem  was  to  identify  the  system  parameters  from  only  input  and  output 
measurements.    The  unknown  covariant  HA  >    tensor  was  chosen  to  be 


[HVi]  - 


.2  -.4      .03  -.7  0.0 

.5  .35      .11  .9  0.0 

.01  1.3  -.33  .7  0.0 

.43  .81  -.05  .4  0.0 

0.0  0.0      0.0  0.0  0.0 


(3.42) 


Therefore,  the  output  equation  consists  of  16  nonzero  terms,  ranging  from 
2x'0)(k)x(°)(k-l)  up  to  .4x(s)(k)x(3'(k-l).  This  model  will  be  called  System  I  in  the 
sequel. 

The  other  type  of  nonlinearity  tested  was  one  that  was  known  to  require 
an  infinite  number  of  terms.  The  particular  example  chosen  for  this  was  the  unit 
step  function  which  simulates  a  saturation  type  of  nonlinearity.  It  is  a  convenient 
choice  since  an  analytical  solution  can  be  calculated. 

The  unit  step  function  has  different  H^  A  tensors  depending  on  the 
order  of  the  model  chosen.  This  is  a  result  of  the  chosen  coordinates  not  being 
orthogonal  (this  will  be  clarified  in  Section  D.)  The  second,  third  and  fourth 
order  (nonlinearity  order)  models  are  given  by 


H 


-    -    0 
2      4 


5    .75    0.0 


H  . 


1.40625    0.0     -1.09375 


.431)! 


H> 


1  _L1 

2      32 


•>•) 


5    1.40625    0.0    -1.09375    0.0 


(3.43c; 
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Notice  that  the  unit  step  function  has  no  memory.  The  expansion  contains  only 
odd  powers  of  x(k)  and  a  constant  term  since  the  unit  step  function  can  be 
written  as  a  constant  plus  an  odd  function  (the  signum  function.)  We  will  refer 
to  these  three  models  collectively  as  System  II  in  the  sequel. 

The  direct  solution  of  the  normal  equations  was  tested  on  these  two 
models  as  was  the  LMS  algorithm.  To  date  the  RLS  algorithm  has  not  been 
verified. 

a.     Direct  Solution  of  Normal  Equations 

To  verify  the  direct  (matrix  inversion)  method  of  solution  of  the 
normal  equations,  a  FORTRAN  program  was  written  which  estimated  the 
correlation  tensors,  and  performed  the  required  matrix  inversion.  The  program 
was  written  so  as  to  allow  the  number  of  points  used  to  approximate  the 
correlations  to  be  varied.  The  maximum  power  of  x(k)  was  also  made  to  be 
adjustable  so  that  the  effect  of  over  or  under  modelling  could  be  studied.  The 
final  adjustable  parameter  was  the  magnitude  of  the  uniform,  zero  mean,  white 
noise  that  was  used  to  excite  the  system.  Adjusting  this  last  parameter  affects  the 
range  over  which  the  resultant  model  is  valid.  In  an  actual  application, 
something  about  the  range  of  expected  system  inputs  would  have  to  be  known  in 
order  to  select  a  good  value  for  this  parameter. 

The  results  for  System  I  are  given  in  Table  3.1  for  several 
combinations  of  the  three  variable  parameters.  The  results  are  remarkably  good 
even  with  as  few  as  30  input  samples.  Since  no  measurement  noise  was 
introduced  this  is  perhaps  not  surprising. 

Overmodelling  did  not  present  any  problems.  The  additional  terms 
were  identified  by  the  algorithm  to  be  essentially  zero.  This  is  evident  in  Table 
3.1c.  Undermodelling  did  introduce  some  inaccuracy.  The  coefficients  identified 
are  significantly  different  from  the  ones  in  the  unknown  system.  However,  the 
coefficients  identified  should  represent  the  best  second  order  approximation  to  the 
system,  which  will  not  in  general  be  the  same  as  the  second  order  coefficients 
contained  in  the  third  order  model. 
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a.  30  data  points  were  used 

The  maximum  power  of  x(k)  was  3 
The  noise  was  uniform  on  (-1,1) 


H 


XgXj 


.200  -.400  .030  -.700 

.500  .350  .110  .900 

.010  1.300  -.330  .700 

430  .810  -.050  .400 


b.  500  data  points  were  used 

The  maximum  power  of  x(k)  was  3 
The  noise  was  uniform  on  (-10,10) 


!H, 


.19993     -.39952 
.49990       .35014 
.010006      1.3000 


.030016      -.70001 

.11001         .90000 

-.33000       .70000 


.43000       .81000      -.050000      .40000 


500  data  points  were  used 

The  maximum  power  of  x(k)  was  4 

The  noise  was  uniform  on  (-1.1) 


[H> 


.20000 

-  .40000 

.03000 

-.70000 

.50000 

.35000 

.11000 

.90000 

1 
0  1' 

.010000 

1.30000 

-.33000 

.70000 

.43000 

.81000 

-.050000 

.40000 

-.48140E- 

06 

-.43483E- 

06 

.34339E-05 

.52571E-06 

26356E- 

06 

.10016E- 

-06 

33811E- 

05 

.83230E 

-07 

.43294E 

-05 

d.  500  data  points  were  used 

The  maximum  power  of  x(k)  was  2 
The  noise  was  uniform  on  (-1,1) 


H 


A  ,,A  , 


.20616 
.75931 
.033444 


-  .80030  .037553 
1.5108  .069939 
1.6619    .24457 


TABLE  3.1:  System  I  Results  Using  Full  Matrix  Inversion. 
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The  results  for  System  II  are  shown  in  Table  3.2.  Significantly  more 
points  were  required  to  accurately  identify  System  II  than  were  needed  for 
System  I.  The  unit  step  function  cannot  be  represented  exactly  by  a  finite 
number  of  polynomials  so  it  is  not  surprising  that  the  solution  is  not  precise. 
There  is  another  factor  that  contributes  to  the  poor  accuracy.  The  choice  of 
functions  used  as  coordinates,  equation  (3.14),  leads  to  ill  conditioned  equations 
(see  Strang  [Ref.  34:  p.  135].)  This  problem  will  be  corrected  in  Section  D. 
b.     Simulation  Using  LMS  Algorithm 

To  verify  the  LMS  algorithm  a  FORTRAN  program  was  written, 
which  used  equation  (3.30),  to  adaptively  identify  the  covariant  tensor  HA  x  . 
The  program  allowed  the  convergence  parameter.  //,  and  the  magnitude  of  the 
uniform,  white  noise  to  be  varied.  The  results  of  the  simulations  are  presented  in 
Table  3.3  and  3.4  for  System  I  and  II  respectively.  Equation  (3.33)  was  used  to 
bound  the  convergence  parameter,  n .  The  input  excitation  noise  was  chosen  to 
be  uniform  on  (-1,1)  resulting  in  a  bound  for  the  convergence  parameter  of 
0</i<0.3558. 

In  general  convergence  was  slow.  The  linear  and  "close  to  linear" 
terms  were  identified  most  rapidly.  The  highly  nonlinear  terms,  involving  high 
powers  of  x(k)  or  x(k-l)  and  their  cross  terms,  were  last  to  be  identified.  Their 
accuracy  never  reached  that  of  the  lower  order  terms.  The  algorithm  was  very 
sensitive  to  the  setting  of  the  convergence  parameter,  p.  A  value  smaller  than 
the  bound  predicted  above  was  used  to  achieve  satisfactory  performance. 

D.    GENERALIZED  COORDINATE  SYSTEMS 

In  Section  B.  a  coordinate  system  was  introduced  that  was  closely  related  to 
the  Volterra  series  (equation  (3.14).)  This  system  was  subsequently  used  in  the 
remainder  of  Section  B  and  in  Section  C.  There-  is  really  little  motivation  for 
choosing  this  particular  set  of  coordinates.  In  fact  there  are  very  compelling 
reasons  to  search  for  other  sets  of  coordinates.  The  set  (3.141  can  lead  to  poorly 
conditioned  sets  of  equations  (see  Strang  [Ref.  34:  p.  135]).  a  fact  that  was 
mentioned  in  the  last  section.  In  higher  order  systems  this  can  become  a  serious 
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a.  1000  data  points  were  used. 

The  maximum  power  of  x(k)  was  3. 
The  input  noise  was  uniform  on  (-1,1) 


!H> 


.49278 

1.4301 

.000519 

-1.1389 

-.037083 

-.10437 

.067337 

.18706 

A0A,1  - 

-.025155 

-.072552 

-  .035273 

.10310 

.084649 

.080036 

-.13334 

-.15251 

b.  15000  data  points  were  used. 

The  maximum  power  of  x(k)  was  3. 
The  input  noise  was  uniform  on  (-1,1). 


H 


A„A,j    - 


.50380 

1.4131 

-.0093646 

-1.1029 

-.0064558 

-.024451 

.030395 

.026762 

-.0016623 

-.021414 

.0044444 

.030779 

00015024 

.024709 

-.029871 

-.01854^ 

c.  500  data  points  were  used. 

The  maximum  power  of  x(k)  was  2. 
The  input  noise  was  uniform  on  (-1,1] 


H, 


,A,      - 


.47416 

.7774" 

.00063628 

0049366 

-.0018654 

0025596 

^020319 

-.10698 

.076991 

d.  5000  data  points  were  used. 

The  maximum  power  of  x(k)  was  2. 
The  input  noise  was  uniform  on  (-1.1). 


U 


.49480         .75558         .0023054 

.017591        .0011427        .024714 

.(25590      .0019598  032117 


TABLE  3.2:  Svstem  II  Results  Using.  Full  Matrix  Inversion. 
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The  maximum  power  of  x(k)  was  3. 

The  noise  was  uniform  on  (-1,1). 

The  convergence  factor,  fi.  was  chosen  to  be  .2. 


a.  After  100  Iterations 


[H, 


b.  After  300  Iterations 


H; 


c.  After  500  Iterations 


[Ha, 


A.J   - 


d.  After  1000  Iterations 


II 


•*[<.■"  i 


e.  After  1800  Iterations 


H.o, 


.21191  -.54901  .23961  -.58084 

.49735  .55747  -.22910  .62263 

.11778  1.3372  -.42404  .85335 

.43573  .41638  -.29008  .45359 


.19133  -.45457 
.52472   .65421 
.032839   1.2483 


-.089402  -.69479 
.055704   .71551 
-.34090   .75055 


.44209   .55546   -.086633   .53219 


.18995  -.35146 

.51438  .55292 

.028302  1.2505 

.43507  .60299 


.20060  -.39481 

.50360  .46156 

.017045  1.2743 

.43879  64136 


.088102  -.68228 

.047074  66728 

.28989  .73386 

082106  .56355 


.028912  -  .70631 

.10901  .74664 

-.32721  75039 

-.070640  .64  572 


.19456  -.38962  .030956  .70562 

.48923  .43032  .11714  .75506 

.0079652  1.2786  -.33206  .7293(1 

43295  .67660  -.053381  .61097 


TABLE  3.3:  System  I  Results  Using  LMS  Algorithm 
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The  Maximum  power  of  x(k)  was  3. 

The  input  noise  was  uniform  on  (-1,1). 

The  convergence  factor,  /*,  was  chosen  to  be  .15. 


a.  After  100  Iterations 


IH, 


27229 

.64437 

-.029974 

-.30266 

.021558 

-.13134 

-.16022 

-.10782 

.062328 

.049563 

-.069966 

-.17809 

26624 

-.00231392 

-.093028 

-.048308 

b.  After  300  Iterations 


(H; 


.49886 

.96021 

-.11685 

-.88415 

.010596 

.10195 

-.060408 

-.12604 

036964 

.15373 

-.17940 

-.32627 

.14847 

.059296 

-.049165 

-.10986 

c.  After  500  Iterations 


H; 


.54880  1.3925  .010187 

-.091216  .14549  .10978 

-.081292  .23020  -.14849 

.0091398  .021484  012741 


-.88497 

-  .0030069 

-.33871 

-.071461 


d.  After  1000  Iterations 


!H; 


.53253  1.5226 

.13698  -  .060273 

.16666  .16941 

068181  .018017 


-.044331  1.0004 

.15947  .036125 

-.096854  -.354  39 

.084221  .0062654 


e.  After  1700  Iterations 


*o  1 


48682  1.3719  -  .058312  .95098 

13368  .025525  .012328  .037707 

.096072  .20474  -.031694  -.21534 

.85853  .047710  -.033043  .025036 


TABLE  3.4:  System  II  Results  Using  LMS  Algorithm 
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problem.  It  is  shown  in  this  section  that  proper  choices  of  coordinates  lead  to 
diagonal  correlation  matrices  so  that  no  costly  matrix  inversions  are  required  to 
solve  the  normal  equations.  This  makes  the  identification  process  almost  trivial. 
It  was  stated  in  Section  B.3  that  the  input  observation  vectors  could  be 
considered  to  describe  a  curve  in  a  hypothetical  (p+l)-dimensional  space.  The 
system  output  is  estimated  based  on  this  curve  (equation  (3.17).)  This  chapter 
generalizes  this  idea  to  include  cases  where  we  have  two  curves,  one  depending  on 
present  and  past  inputs,  the  other  depending  on  past  outputs.  Finally,  the 
chapter  concludes  with  several  non-trivial  numerical  simulations. 
1.     Choices  of  Coordinate  Systems 

A  point  in  (p-(-l)-dimensional  space  is  defined  by  the  variables 
x°.  x2,  ...,  xp.  In  Section  B.3  (equation  (3.14))  these  variables  were  chosen  to  be 
parametric  functions  of  the  variable  x(k),  the  system  input.  The  particular 
choice  presented  there  was  chosen  to  correspond  to  the  Yolterra  series.  It  is 
desirable  to  pick  a  system  of  coordinates  that  ensure  a  diagonal  correlation 
tensor,  as  this  allows  solution  of  the  normal  equations  without  matrix  inversion. 
A  tensor  T  is  diagonal  if  its  components  obey  the  rule 


rp-*0 


K>      'or  ^0-  ^  N-  1  •    ^l-^"1      .  •         •   ^  N^l      .-  ^ 


-2~;]  ~_1  (3.44) 

otherwise 


Note  that  only  tensors  of  an  even  order  (possessing  an  even  number  of  indices) 
can  possibly  be  diagonal.  In  general  components  satisfying  the  upper  condition  of 
equation  (3.44)  are  called  diagonal  elements,  or  diagonal  components. 
Components  that  are  not  diagonal  are  called  off  diagonal. 

Two   conditions   are   required   in   order  to  ensure  a   diagonal   correlation 
tensor.    The  first  is  a  result  of  the  following  theorem. 

a.      Theorem  3.1 

If  a  set  of  functions  {  f0(x),  f,(x) 1n(x)}  »re  defined  with  the  property 

that 
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-  (X  V 

where  w(x)  is  a  positive  weighting  function,  then  the  set  of  random  variables.  ZA. 
defined  as 

Z°  =  f0(X) 
Z1  =  ft(X) 

(3.46) 

Zp  =  fppc) 

where  X  is  a  random  variable  with  probability  distribution 

Px(x)  =  Cw(x)  (3.47) 

are  uneorrelated.  In  equation  (3.47)  the  constant.  C.  and  weighting  function. 
w(x).  must  be  chosen  so  that  px(x)  satisfies  the  definition  of  a  probability  density 
function. 

b.     Proof 


«ZAZM|   =    EJf; 


(X)f,(X)  (3.48a) 


=  JfA(x)fA(x)Px(x)dx  (3.48b) 

—  00 

--   JfA(x)fM(x)Cw(x)dj(  (3.48c) 

—  oc 

DO 

-  CjffA(x)f„(x)w(x)dx  (3.48(1) 

—  00 

Substituting  equation  (3.45).  yields 

e/zaz4      (K;,  iorX=  "  (3.49) 

Choosing  a  set  of  functions  in  accordance  with  equation   (3.45)   and 
ensuring  that  f0(x)   is  a  constant  function   (ie:  is  a  constant)   is  the  first  condition 
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that  must  be  satisfied  in  order  to  obtain  a  diagonal  correlation  tensor.  The 
system  must  be  excited  with  samples  drawn  from  a  random  process  with 
probability  distribution  given  by  (3.47).  In  addition,  if  the  random  process  is 
chosen  to  be  zero  mean  and  strictly  white  (ie.  independent  which  implies 
uncorrelated),  then  the  correlation  tensor  will  be  diagonal.  This  last  condition  is 
equivalent  to  requiring  that 

E|x(m)x(n)|  =  <5(m-n)  (3.50a) 

p(x(n),x(m))  =  p(x(n))p(x(m))  (3.50b) 

and 

E{x(m)i  =  0  (3.50c) 


where  x(m)  and  x(n)  are  two  input  samples  taken  at  times  m  and  n  respectively. 
It  is  a  straight  forward  matter  to  show  that  if  these  conditions  are  met.  the 
correlation  tensor  will  be  diagonal. 

The  conditions  presented  above  imply  that  different  sets  of 
coordinate  functions  should  be  used  depending  on  the  probability  distribution  of 
the  noise  used  to  excite  the  system.  Different  noise  distributions  have  the  effect 
of  weighting  the  error  differently.  Consider  that  a  Gaussian  noise  will  contain 
samples  of  all  amplitudes  while  uniform  noise  is  bounded.  If.  for  example,  the 
system  contains  a  saturation  type  of  nonlinearity.  the  uniform  noise  may  not 
detect  its  presence  if  its  maximum  amplitude  i-  not  sufficiently  large. 
Theoretically.  Gaussian  noise  contains  sample>  of  all  amplitudes  and  will  excite 
all  modes.  On  the  other  hand  it  may  be  known  that  rh<  system  input  never 
fxc<>r-d>  a  certain  maximum  value  and  so  h  bounded  input  will  'he  >uitable. 

If  two  models  of  a  system  are  constructed,  in  two  different  coordinate 
systems  but  using  the  same  input  noise  (only  one  >et  can  possibly  lead  to  a 
diagonal  correlation  tensor)  then  the  two  solutions  will  be  equivalent.  One 
solution  can  be  transformed  into  the  other  by  performing  a  change  of  coordinates 
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in  accordance  with  equation  (2.22).  Therefore,  it  always  makes  sense  to  identify 
the  system  using  the  coordinate  system  which  leads  to  a  diagonal  tensor.  Other 
representations  can  then  be  calculated  as  desired.  Solutions  obtained  by  using 
different  input  noise  density  functions  are  not  equivalent  even  if  the  coordinate 
systems  used  were  the  same.  The  effect  of  the  different  distributions  is  to  weight 
the  errors  differently.  A  uniform  input  noise  will  weight  the  errors  equally,  while 
a  Gaussian  noise  will  emphasize  the  importance  of  errors  made  for  small  inputs. 
In  general,  transformations  between  solutions  obtained  using  different  excitation 
noise  distributions  cannot  be  found.  Choosing  an  appropriate  distribution 
requires  knowledge  of  the  expected  system  input  signals. 

The  Hermite  polynomials  lead  to  a  diagonal  correlation  tensor  if  the 
input  is  white,  zero  mean,  Gaussian  noise.  Similarly,  Legendre  polynomials 
should  be  used  in  the  case  of  uniform  noise.  It  is  convenient  to  normalize  the 
coordinate  functions  so  that  the  diagonal  components  of  the  correlation  tensor  are 
all  ones.  To  identify  the  system  parameters,  only  the  cross-correlations  of  the 
right-hand-side  need  to  be  calculated.  This  fact  was  first  discovered  by  Lee  and 
Schetzen  |Ref.  10]. 

2.      Recursive  Models 

Recursive  models  have  been  used  very  successfully  for  modelling  linear 
systems  [Ref.  35].  Among  their  advantages  is  infinite  memory,  and  the  ability  to 
model  a  system  without  knowledge  of  its  input.  The  latter  property  allows  these 
models  to  be  employed  in  such  areas  as  speech  processing  where  input  signal  arc 
difficult  or  impossible  to  measure.  The  assumption  is  made  that  The  input  i^ 
white  noise. 

Recursive  discrete-time  nonlinear  expansions  have  also  been  proposed 
(see  for  example  'Ref.  14.15.16]).  Recursive  models  po-ses  infinite  memory,  and  so 
may  require  fewer  terms  to  accurately  represent  long  memory  systems.  However, 
nonlinear  recursive  models  also  posses  infinite  uonlinearity.  To  understand  why 
this  is  true,  consider  the  following  example. 
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a.     Example  3.2 

Consider  the  following  recursive,  discrete,  nonlinear  system 

y(k)  =   ay2(k-l)  -  u(k)  (3.51) 

where  u(k)  =  b<5(k)  (where  £(k)  is  a  unit  sample)  and  where  a  and  b  are  arbitrary 
constants. 

The  output  (y(k))  of  the  system  for  k  =  1,  ...,  K  is 

y(0)  =  b  (3.52a) 

y(l)  =  ab2  (3.52b) 

y(2)  =   a(ab2)2  -  aV  (3.52c) 


y(K)  =   a2""1  b2*  (3.52d) 

It  is  clear  that  the  nonlinearity  of  the  system  increases  with  time.  Unlike  a  linear 
system,  stability  in  a  recursive  nonlinear  system  is  determined  not  only  by  the 
system  parameter,  a.  but  also  by  the  input  function.  It  is  also  difficult,  in 
general,  to  predict  for  what  classes  of  input  a  particular  system  will  be  stable. 
By  analogy  to  the  linear  case  we  will  refer  to  these  types  of  models  as  Auto- 
Regressive  or  AR. 

It  is  possible  to  also  expand  the  output  as  a  combination  of  both  past 
and  present  inputs  and  past  outputs.  This  type  of  model  was  proposed  by  Parker 
and  Perry  'Ref.  13].  We  will  refer  to  this  type  of  model  as  an  Auto-Regressive. 
Moving- Average  or  ARMA  model.  It  will  have  the  same  stability  problems  as 
does  the  AR  model. 

Equation  (3.17)  can  be  used  to  model  an  AR  nonlinear  system-  if  the 
proper  choice   is  made  for  the  observation  vectors.     Using  the  same  coordinate 

functions  as  given  in  equation  (3.14)  but  using  y(k  -  i).  i  =  1 N,    as  an  input 

parameter,  yields  appropriate  observation  vectors.  That  is. 
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yAi(k-L)  =,  !y(k-i)]W  (3.53) 

defines    the    components   of  the   i-th   observation   vector.     The   model   equation, 
equivalent  to  equation  (3.17),  becomes 

y(k)  =  yAl(k-l)...yAN(k-N)HAi.   An  (3.54) 

where  HA  ..  A     is  again  a  (p+l)-th  order  covariant  tensor  containing  the  system 

parameters.    This  model  is  an  extension  of  the  familiar  linear  autoregressive,  or 
AR  model. 

The  normal  equations  derived  for  the  nonrecursive  case  can  be  used 
to  solve  for  the  HA        A    tensor  in  this  case  as  well.  The  identification  process  is 

described  in  Figure  3.5.    The  assumption  is  made  that  the  system  is  recursive  and 
driven  by  a  white  noise,  u(k).    Its  output  can  then  be  described  by 

y(k)  =  yAl(k-l)...yAN(k-N)Hv     ^  +  u(k)  (3.55) 

This  output  signal  is  delayed  and  fed  into  the  analysis  model.  The  error  signal 
e(k)  is  given  by 

e(k)  =  y(k)  -  y(k)  (3.56a) 


=  y  »(k-l)...y"(k-N)Hv     AN-u(k)-y'(k    ])...v  N(k-N)HAi     .x 


(3.56b) 


When  the  model  parameters  exactly  equal  the  actual  system 
parameters  the  error  signal  will  equal  the  input  white  noise.  For  this  reason  the 
analysis  model  is  often  called  a  whitening  or  bleaching  filter.  The  analysis 
model  is  nonrecursive.  It  uses  past  values  of  the  system  output  to  make  a 
prediction  of  the  present  system  output.  The  normal  equations  (3.25c)  apply  to 
this  situation  with  the  observation  vectors,  x(k-i).  replaced  by  the  vectors  defined 
in  equation  (3.53).  The  normal  equations  for  the  recursive  nonlinear  model  can 
be  written  as 
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Figure  3.5:  Recursive  Model  Identification 
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Ejy  l(k-l)...y  N(k-N)y"l(k-l)-yPN(k-N)|HA]       Xy 

=  E|y(k)y'Jl(k-l)...y'iN(k-N)|  (3-57) 

These   equations    are   a   tensor   equivalent    of  the    Yule-Walker   equations.     The 
corresponding  mean-square  error  is 

E|e2(k)|  =  E|y2(k)|  -  E|y(k)yAl(k- 1)  ■  ■  •  yAN(k-N)JHv      x„  (3.58) 

In  this  case  it  is  not  obvious  that  one  set  of  coordinates  will  yield 
better  results  than  another.  The  correlation  tensor  appearing  on  the  left-hand- 
side  of  equation  (3.57)  cannot  be  guaranteed  to  be  diagonal  since  the  probability 
distribution  of  y(k)  cannot  be  controlled.  Techniques  must  be  devised  which 
choose  the  coordinate  system  "on  line"  as  the  statistics  of  y(k)  are  determined. 
For  example,  in  the  linear  lattice  the  coordinate  system  is  chosen  by 
orthogonalizing  the  sequence  y(k)  using  a  Gram-Schmidt  procedure.  In  this  way 
the  coordinate  system  is  not  determined  until  run  time.  Extension  of  these  ideas 
is  left  until  Chapter  4. 

The  tensor  model  presented  is  also'  suitable  to  represent  a  nonlinear 
ARMA  model.  This  type  of  model  takes  into  account  all  available  information 
(input  and  output)  and  so  should  be  more  accurate.  It  may  also  lead,  in  some 
cases,  to  a  lower  order  solution  than  either  the  AR  or  MA  model.  An  ARMA 
tensor  model  can  be  written  as 

y(k)  =  xAu(k)...x"M(k-M)yM,(k-l)...vMN(k-N)HAo.     ,mMi       „„  (3.59) 

Relations   analogous  to   (3.57)    and   (3.58)    for   the  normal   equations 
and  the  minimum  mean-square  error    for  the  ARMA  model  can  be  obtained  in  a 
straight  forward  manner. 
3.      Simulation  Results 

The  concepts  developed  in  Section  D  of  this  chapter  were  verified  using 
FORTRAN    simulations.      The    coordinate    functions    were    chosen    to    be    the 
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normalized  Legendre  polynomials  because  of  the  ease  in  generating  uniform  noise 
on  a  digital  computer.  By  assuming  that  the  correlation  matrix  was  diagonal, 
the  solution  to  the  normal  equations  was  determined  without  matrix  inversion. 
This  solution  was  verified  by  generating  the  correlation  on  the  left-hand-side  of 
the  normal  equations  and  performing  the  required  matrix  inversion.  The  LMS 
adaptive  algorithm  was  also  tested  using  the  Legendre  polynomials. 
The  Legendre  polynomials  used  are  given  by 


fo(x)  =   1.0 


(3.60a) 


f,(x)  -v^x 


(3.60b) 


f2(x) 


x'2)  -  ±\ 


(3.60c) 


fi(x) 


175 


;3.60d) 


These  functions  are  normalized  so  that  the  correlation  tensor  will  have  ones  along 
the  diagonal  when  excited  with  zero  mean  white  noise  which  is  uniform  on  (-1,1). 
The  simulation  models  used  were  similar  to  those  used  in  Section  C.  The 
coefficients  for  System  I  were  identical  to  those  given  in  (3.42).  although  thi- 
time  they  are  coefficients  of  Legendre  polynomials  so  they  do  not  represent  the 
same  system.  System  II  was  again  a  unit  step  function  which  has  a  tensor 
representation,  in  terms  of  the  Legendre  polynomials  given  by  [Ref.  26:  p.  1526], 


HA, 


4 


\  175 

80 


.61a' 


0.5    0.433013    0.0       0.165359 


(3.61b 


The  results  of  the  simulations  for  System  I  are  presented  in  Tables  3.5 
and  3.G.  Results  for  System  II  are  given  in  Tables  3.7  and  3.8.  In  all  cases  it  is 
obvious  that  the  simulation  results  are  very  good. 
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15000  sample  data  points  were  used. 
The  maximum  power  of  x(k)  was  3. 


a.  Solution  Without  Using  Matrix  Inversion 


iH 


V,J 


.19158 

-.38475 

.0041187 

-.67056 

.49067 

.38457 

.089293 

.93948 

.012253 

1.2925 

-.35452 

.71923 

.41049 

.85202 

-.074247 

.43011 

b.  Solution  Using  Matrix  Inversion 


H 


Vi 


.20000 

-.40000 

.030000 

-.70000 

.50000 

.35000 

.11000 

.90000 

.01000 

1.3000 

-.33000 

.70000 

.4800 

.81000 

-.050000 

.40000 

TABLE  3.5:  System  I  Model  Parameters  Obtained 

a.  Using  No  Matrix  Inversion,  and 

b.  Using  Matrix  Inversion. 
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The  maximum  power  of  x(k)  was  3. 

The  convergence  factor  was  chosen  to  be  .01. 


a.  After  200  Iterations 


H 


A0Aj 


.14678   -.45852 
.42720    .25876 
-.077103   1.1898 


-.036398  -.72277 
.018320   .83079 
-.45388   .61647 


.34264 


.68052   -.16955   .30546 


b.  After  500  Iterations 


!H, 


.20013  -.40047  .030509  -.70066 

.50038  .35077    .11136    .90025 

.0099732  1.2996  -.32958   .70026 

.43008  .81014  -.050063   .39875 


c.  After  1000  Iterations 


H 


A„A,  - 


.20000  -.40000  .030007  -.70000 

.50001  .35000  .11001  .8999 

.010008  1  3000  32999  .69999 

.43001  .81000  -.049986  .40000 


TABLE  3.6:  System  I  Model  Parameters  Using  LMS  Algorithma. 
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a.  500  data  points  were  used. 

The  maximum  power  of  x(k)  was  3. 


No  Matrix  Inversion  Solution 


H 


AgA, 


.52104 
-.011098 
.0011160 

.040052 


.43930  .0090070  -.17515 

-.0046220  .014464  .020122 

-.010606  -.010466  .021310 

-.0071009  .058263  .052566 


Using  Matrix  Inversion 


'H, 


.50157 

-.0084508 
.0017041 


.43054 
.0032387 
-.0009167 


-.0066652    -.0024676 


-.011153 

-  .0000970 

.0062532 

.012192 


-.16530 
-.0006363 
-.0059126 
.0014058 


b.  15000  data  point  were  used. 

The  maximum  power  of  x(k)  was  3. 


No  Matrix  Inversion  Solution 


!H, 


.50050 
.0015150 
-.0071103 
.0022021 


.43363 

.0003272 

-.0004279 

.0018209 


.0043279 
-.0012861 

.012933 
-.0021332 


-.16772 

.0003738 

014358 

.0060674 

Usins  Matrix  Inversion 


|Ha„a. 


..".0089 
00)8213 
0036401 

0003990 


-.0005741 
.0008143 
.0001746 


-.0001744 
-  .0014410 
.0043811 
- .0003974 


.16722 

.0000789 

-.0008015 

.0000407 


o.t.  System  II  Model  Parameters  Obtained 
a.  Using  No  Matrix  Inversion,  and 
1).  Using  Matrix  Inversion. 
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The  maximum  power  of  x(k)  was  3. 

The  convergence  factor,  /*,  was  chosen  to  be  .005. 


a.  After  300  Iterations 


IH, 


.46464 

.011790 

-.010943 

.0014149 


.41356  .0066999      -.16385 

.012146  .0067764  -.0093346 
.014180       .020212        .018152 

.0065302  .0093879  -.0047952 


b.  After  500  Iterations 


AnA,; 


.50556 

.43020 

-.00098824 

-.16796 

.015344 

-.010684 

.020836 

.0098449 

.0082813 

-.0052115 

.0090727 

-.010251 

012877 

-.0094530 

-.0046525 

.0072316 

c.  After  1000  Iterations 


.49868    .44629  .0053940  -.17370 

.020577  .0091397  -.010215  .0065185 

.012857  .020072  -.0064996  -.010541 

.0020197  .0031219  -.0064275  -.0037119 


c.  After  1800  Iterations 


H,  , 

'  0'  1 


.50115     .43727 

.0000129  -.0063172 

0065666  .0084999 

-0075634  .0058037 


.0059911  -.16508 

.012713  .010275 

.0068103  -  .010135 

-.15156  .018235 


TABLE  3.8:  System  II  Model  Parameters  Using  LMS  Algorithm 
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IV.    REVIEW  OF  LATTICE  FILTER  STRUCTURES 

This  chapter  reviews  lattice  filter  theory.  These  filters  were  first  proposed  in 
connection  with  the  linear  prediction  of  speech  by  Itakura  and  Saito  [Ref.  36]  in 
1971.  They  developed  a  new  approach  based  on  a  partial  correlation  (PARCOR) 
coefficient.  Since  that  time  the  filter  structure  they  proposed  has  come  to  be 
known  as  a  Lattice  (or  sometimes  also  called  Ladder)  filter.  The  properties  of  the 
filter  have  been  exhaustively  studied  by  many  researchers  [Ref.  37].  Lattice 
filters  have  been  successfully  applied  to  problems  in  various  disciplines. 

The  great  interest  in  the  lattice  approach  stems  from  it's  property  of 
orthogonality.  This  property  allows  the  filter  to  be  updated  in  order,  without 
recalculation  of  all  the  previous,  lower  order,  filter  coefficients.  Orthogonality  also 
leads  to  a  nice  physical  structure,  a  cascade  of  first  order  sections,  and  so  is 
appropriate  for  efficient  hardware  or  software  implementations.  Finally,  it  has 
been  shown  that  the  lattice  owes  its  robust  numerical  behaviour  to  this  property 
of  orthogonality  [Ref.  38:  pp.  128-136]. 

This  chapter  begins  with  a  derivation  of  the  one-dimensional  (1-D)  lattice. 
The  approach  taken  in  the  derivation  is  somewhat  novel  in  that  it  begins  by 
expressing  the  linear  prediction  in  terms  of  an  uncorrected  error  sequence.  This 
is  regarded  a^  a  change  of  coordinate  systems  and  used  to  develop  the  Levinson 
algorithm  and  the  lattice  filter.  In  the  following  section  generalized  forms  of  the 
Levinson  algorithm  and  lattice  are  derived.  It  is  shown  that  these  also 
correspond  to  coordinate  transformation.  Next,  the  Schur  algorithm,  which  is  a 
method  for  generating  the  required  filter  coefficients  (Lattice  parameters) 
directly,  given  only  a  knowledge  of  the  correlation  matrix,  is  reviewed.  In  the 
following  chapter  the  generalized  lattice  filter  is  used  to  develop  new. 
multidimensional  (specifically  2-D)  lattice  filters.  In  Chapter  G.  a  new  nonlinear 
lattice. filter,  based  upon  this  generalized  lattice  formulation,  is  presented. 
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A.    1-D  LINEAR  AUTOREGRESSIVE  LATTICE  FILTER 

We  define  the  N-th  order  forward  error  sequence  of  an  autoregressive  model 
as 


eN(k)  =  y(k)  -    £  hAN  y(k  -  A)  (4.1) 

A=l 

This  can  equivalently  be  written  in  tensor  notation  as, 

eN(k)  =  hAV  (4.2) 

where 

[yA]T=  |y(k)    y(k-l)    ■■■    y(k-K)]  (4.3) 

and 

[hAN|  =  [1    -h,N    -h2N    ■••    -hNN0...0]  (4.4a) 

where  for  convenience,  we  have  made  all  vectors  of  length  K  +  l  (ie;  A  =  0.  ....  K). 
where  K  is  some  arbitrary  maximum  length.    Note  that, 

h0N  =  1  "  (4.4b) 

The  y;  can  be  considered  to  comprise  a  coordinate  system  in  an  K+l 
dimensional  space.  The  vector  ,yA;  represents  a  single  realization  from  a  random 
process.  As  mentioned  in  the  Chapter  I,  there  will  in  general  be  many  such 
vectors  corresponding  to  all  the  possible  realizations  of  the  random  process. 

Because  the  error.  eN(k).  is  a  scalar  (invariant)  it  must  remain  unaltered 
regardless  of  the  choice  of  coordinate  system.    We  may.  therefore,  write 

-N(k)  =   KAV  (4.5) 

where  the  vA    are  defined  as 
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-A': 


y(k) 

y(k-l) 
y(k-2) 
y(k-S) 


hoV(k-l) 
h12y(k-2)h02y(k-i; 


y(k-K)  -hKK_-2V(k-K+l) 


hoK_,y(k-i) 


(4.6) 


and 


\K?A 


,y]-  [1    -Kj    -K2   ••■    -KN0...0]  (4.7) 

where  the  Vs  are  chosen  so  that  the  components,  yA  .  are  uncorrelated.  The 
stochastic  form  of  the  Gram-Schmidt  procedure  can  be  used.  A  discussion  of  this 
method  can  be  found  in  [Ref.  39:  pp.  382-383].    By  uncorrelated  we  mean 


EV> 


I  ^  a  •    for  A  '  =  fi 

\0        otherwise 


(4.8) 


The  reader  familiar  with  lattice  structures  will  recognize  the  components  of  yA  as 
the  backward  prediction  errors.  That  is  they  are  the  errors  in  predicting  y(k-N) 
from  the  next  N-l  values:  y(k-N-f  1).  ....  y(k-l). 

It  is  a  straight  forward  matter  to  solve  the  prediction  problem  given  by 
equation  (4.5)  (ie;  solve  for  the  K's)  because  of  the  uncorrelatedness  of  the  chosen 
coordinate  system.  Using  an  approach  similar  to  that  presented  in  Chapter  III. 
Section  C.  a  set  of  normal  equations  can  be  formulated  for  this  problem.  In  this 
case,  however,  the  correlation  matrix  is  diagonal  so  that  there  is  no  inversion 
necessary  to  obtain  the  solution.  Minimizing  the  mean  square  value  of  the  error. 
eN(k).  given  by  (4.5).  with  respect  to  the  K*-"s.  yields 


K,N 


K ■v(k)>'    • 

Ef',m) 


(  ...   S-) 

y(k)\     > 


4.9] 


Having  obtained  a  solution  to  the  orthogonal  problem,  we  wonder  if  it  cannot 
be    employed    to    advantage    to    simplify    the    calculations    required    to    solve    the 
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original    autoregressive    problem.      From    pquation    (2.13).    we    know    that    the 
relationship  between  the  old  and  new  components  can  be  expressed  as 


9^-v 


9yA 


(4.10) 


where 


ay 


X' 


3y; 


0 

1 


K-2 


0 
0 

1 

h,2 


K-2 


0     _hK-l     _"hK-l 


0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

K-l 
lK-2 

1 

4.11 


Now.  since 

eN(k)  =   yV  =   yvKAN 


yAllL    K,N- 


9y' 


(4.12a) 


(4.12b) 


th 


en 


i  n      r.-  n  9y 

dv 


4.i: 


Equation  (4.13)  gives  the  relationship  between  the  autoregressive  model 
parameters  and  the  KAN-'s.  This  result  will  be  used  in  the  proof  of  LevinsorTs 
algorithm. 

1.      Levinson-Durbin  Algorithm 

In  1947  in  his  classic  paper  [Ref.  -JOl  Levirison  developed  a  method  for 
recursively  solving  the  normal  equations.  Beginning  with  a  zero  order  solution 
successively  higher  order  solutions  are  calculated.  This  algorithm  can  be  used  to 
exploit  the  Toeplitz  nature  of  correlation  matrices  of  stationary  random  processes 
in  order  to  reduce  the  required  number  of  computations.  The  algorithm  as 
presented  in  this  work  is  actually  a  simplified  version  of  Levinson's  original  from. 
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It  assumes  that  the  equations  being  solved  are  the  Yule- Walker  equations  ([Ref. 
41],  see  also  Chapter  3.  Section  D)  so  that  the  right-hand-side  of  the  equations 
contains  only  terms  that  will  be  present  in  the  correlation  matrix  (left-hand-side) 
of  the  next  higher  order  model.  This  simplification  is  due  to  Durbin  [Ref.  42]. 
We  will  refer  to  this  algorithm  often  simply  as  the  Levinson  algorithm,  however, 
Durbin 's  contribution  is  acknowledged. 

a.  Theorem  4.1  Levinson's  Algorithm  (Durbin's  Form) 

The  autoregressive  model  parameters  may  be  calculated  recursively 
from  the  relation 

ihr1]  =  [hAN]  -  KN.,S!hAN]  (4.14a) 

where 

!hANj  =  !-h*    -h,N    •••    -hNN_,    10   •■■   0]  (4.14b) 

and 

SfhAN]=  10    -h0N    -h,N    •••    -h£_,    1   0   ••   0]  (4.14c) 

The  operator  S  has  the  effect  of  shifting  the  components,  h/  one  position  To  the 
right.  Note  that 

For  stationary  processes  the  following  simplification  applies 

h/1  =  *»n-a-]      for  A  -  0 N  (4.14d) 

b.  Proof  (by  Induction) 

Using  equation  (4.13)  for  the  first  order  model  we  have 

hi    =     ho1    h,1    0    •  •  •    0  (4.15a) 

1       K,    0    ■  ■  •    0  (4.15b) 

1    0    ■  •  ■    0        K,  0    10-0  (4.15c) 

=  [hA°j  -  K.sfh,0  (4.15d) 
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so  that  equation  (4.14a)   holds  for  N=l.     We  assume  the  recursion  holds  for  the 
N-th  order  model 


[hANl  =   ih^-1]  -  KNSjh 


N-l 


From  equation  (4.13)  we  have 

A' 

—      A,  A'  =  0,  ...,K 


h>    -  Kv  — - 

3y 


Therefore, 


,hA       - 


-K,  +  K2h0]  +  K3h02  +...+  KNh 


N-l 


K2  +  K3h;  +  K4h,  +...+  KNhi 


N-2 


KNhNJ2 


Now.  for  the  (N  +  l)-st  order  model 


hi 


.N-l 


IT 


1 


K,  -  K2h0'  -  K,h0"  +...+  KN^h0* 
K2  -  K3h,a  +  K4h"  +...+  KN+1hrN 


KN  -  KN,,hN_  , 


(4.16) 


(4.17) 


(4.18) 


'4.19a 
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■K,  -  Kji,,1  +  K8hc 


Kvh 


N-l 


N"0 


■K, 


K3hf  +  K4h,3  +  ...+  KNh,J 


N-l 


KN 


—    KN+] 


hn* 


V.  N 

1 

0 


(4.19b) 


hAN    -  KN+1S!hAN] 


(4.19c) 


And  so  the  desired  result  has  been  confirmed. 

The  proof  presented  does  not  require  that  the  process  be  stationary 
and  so  the  condition  of  (4.14c)  is  not  necessary.  One  is  then  faced  with  the 
problem  of  how  to  determine  the  hAN's.  This  question  is  answered  in  the  next 
section  where  a  generalized  form  of  the  algorithm  is  derived.  We  note  that  the 
condition  of  (4.14c)  implies  the  following 

(4.20) 


h^" 


KN  =  hc 


or  that  the  second  column  of  the  matrix  given  in  equation  (4.11)  contains  all  the 
lattice  coefficients. 
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The  following  significant  observation  about  equation  (4.14a)  is  made. 
The  transpose  of  the  term  within  the  second  bracket  on  the  right  hand  side 
corresponds  exactly  with  the  rows  of  the  coordinate  transformation  matrix  of 
equation  (4.11).  The  Levinson  algorithm,  therefore,  represents  a  recursive 
method  for  finding  an  orthogonalizing  coordinate  transformation. 

2.     The  1-D  Lattice  Structure 

It  has  been  shown,  in  the  previous  section,  that  the  autoregressive  model 
parameters  can  be  calculated  in  a  recursive  manner  using  the  Levinson  algorithm. 
It  was  also  shown  that  a  model  could  be  built  in  an  orthogonal  coordinate  system 
where  the  model  parameters  were  given  by  the  KA-'s.  The  question  arises,  can  a 
filter  structure  be  devised  which  represents  y(k)  in  the  orthogonal  coordinate 
system?  The  answer  is  affirmative^  It  will  be  shown  in  this  section  that  the 
Lattice  filter  is  the  required  structure  for  the  stationary  case.  A  more  general 
solution  is  presented  in  Section  B  of  this  chapter. 

The  desired  result  is  obtained  in  a  straight  forward  manner  by 
multiplying  both  sides  of  equation  (4.14a)  by  yA.  This  yields 

eN-!(k)  =  eN(k)  -  KNV1rN(k-N-l)  (4.21) 

where  the  quantity  rN(k-N-l)  =  yA  evaluated  at  A'  -  N+l.  As  was  mentioned 
earlier  in  this  chapter,  the  quantity  rN(k-X-l).  is  generally  known  as  the 
backward  prediction  error  since  it  corresponds  to  the  prediction  of  the  point  y(k- 

N-l)  from  the  N  future  points  y(k-l) y(k-N). 

Assuming  stationarity.  we  can  use  (4.14c)  and  (4.14a)  to  obtain 

fhAN+1i  -  SfhAN]  -  KN.rhAN'  (4.22) 

which  leads  directly  to  the  equation 

rN+1(k-N-l)  =  rN(k-N-  1)  -  KN  +  1eN(k)  (4.23) 

Equations  (4.21)  and  (4.23)  define  the  Lattice  form  of  the  whitening  filter  (see 
Chapter  3,  Section  D).  This  is  also  sometimes  referred  to  as  the  analysis  model. 
The  structure  is  illustrated  in  Figure  4.1. 
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In  order  to  develop  a  synthesis  model  we  need  only  to  solve  equation 
(4.21)  for  eN(k)  since  we  know  that  e°(k)  is  equivalent  to  v(k).  Therefore,  the 
synthesis  equations  are 

eN(k)  =  eN+1(k)  +  KN+1rN(k-N-l)  (4.24a) 

rN+1(k-N-l)  =  rN(k-N-l)  -  KN+1eN(k)  (4.24b) 

The  resulting  structure  is  as  illustrated  in  Figure  4.2.  Note  that  in  this  model, 
y(k)  is  indeed  being  expressed  as  a  weighted  sum  of  the  backwards  errors,  which 
represent  an  orthogonal  coordinate  system. 

This  concludes  the  discussion  of  the  classical  Lattice  formulation. 

B.  GENERALIZED  ORDER  UPDATE  RECURSIONS 

In  this  section  a  more  general  linear  prediction  problem  is  considered.  No 
assumptions  are  made  as  to  the  origin  of  the  data.  In  fact,  the  data  need  not 
represent  a  time  series  at  all,  and  certainly  shift  invariance  is  not  required.  The 
ordering  of  the  data  is  simply  chosen  in  some  convenient  fashion.  The  generality 
of  this  formulation  will  allow  its  application  to  multidimensional  and  nonlinear 
problems. 

1.      Definitions  and  Formulation 

In  this  section  quantities  are  defined  that  will  be  required  to  complete 
the  statements  and  proofs  contained  in  the  remainder  of  the  chapter. 

A  realization  of  the  random  process  Y  is  given  by  the  vector 
jyAj,  0  <;  A  v  K. 

The  error.  ekN_,.  in  predicting  the  element  yk"'  from  the  previous  N 
elements  of  Y  is  given  by 

ekN;,  =   hAN(k+lV       for  A  -   0 K  (4.25) 

where 

hAN(k+l)j  =  [0   ■••   0     -h£N+1    -hkN_N+2    ■■■    -hkN    1    0   ••■    0]  (4.26) 
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r  3(k-3) 


r°(k)  r1  (k-1)  r2(k-2) 


Figure  4.1:  1-D  Lattice  Analysis  Model. 


*r-*-e3(k) 


r  3(k-3) 


r°  (k)  r1  (k-1)  r2(k-2) 


Figure  4.2:  1-D  Lattice  Synthesis  Model. 
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The  norm  of  this  prediction  error  is  given  by 


;k  +  l  i 


N     \2l 


Wrt) 


1/  2 


A  normalized  version  of  the  forward  error  is  given  by 


«N    - 


ek  +  l 


ek+i 


(4.27) 


(4.28a) 


aAN(k+l)y; 


(4.28b) 


where 


aAN(k+l) 


_  N        I       ! 

ek  +  l 


h,N(k-l) 


(4.29) 


The   backwards   prediction   error,  rk1N,   is   the   error   associated  with   the 
prediction  of  yk_N  given  the  next  N  elements  of  Y.    It  is  given  by, 

rkN-N  =  hAN(k-NV  (4-30) 

where 

fhf(k-N)]  =  [0   •  •  •   0   1    -h£N+J    -hkNN+2    •  •  •    -hkN   0   •  •  ■   0]  (4.31) 

Again,  a  normalized  version  can  be  defined  as 


_N  Tk-N 

rk-N  - £ — 

rk-  N 


(4.32a) 


-   bAN(k-N)yA 


;4.32b) 


w 


here 


b? 


hr(k-N) 

N 
rk_N  j 


and 


rk-N 


B[(rkN-N)2 


\l  2 


(4.33) 


(4.34) 
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In  order  to  generalize  the  results  of  Section  A  of  this  chapter  we  need  to 
introduce  two  different  families  of  coordinate  systems.  We  will  refer  to  members 
of  these  families  as  either  forward  or  backward  local  coordinates.  The  term  local 
is  used  because  a  different  coordinate  system  will  be  associated  with  every  value 
of  k  and  N.  We  define  the  (k  +  1)  -  indexed,  N-th  order  forward  coordinate 
system  as 


[yA'(k+l,N)j  = 


,k-N 


,k-N  +  l 


N-l         k-N-r2 


-  K-S^V 


N-l  ,k 


hk     y 


,k-l 


k  +  1 


l,,k 


hkV 


(4.35) 


The    corresponding    coordinate    transformation    from    the    unprimed    coordinate 
system  to  this  local  forward  coordinate  system  is  given  by 

By-1  (k-l.N) 
BvA 


88 


1 

-hk-N- 

N-l 
•2 

I  N-1 
nk-N  +  3 

"  hk-l 

v.  N   ' 

0 

1 

7   N-2 
_nk-N  +  3 

-hk-i 

-hk 

0 

0 

1 

r  n-s 
-hk-i 

"hk 

-L1 


0 


0 


I 


(4.36) 


where  the  o's  are  zero  matrices  and  the  I's  represent  identity  matrices. 

Similarly,  we  define  the  (k-N)-indexed,  N-th  order  backward  coordinate 
system  as 


.A" 


yA    (k-N.N)   = 


k-K 

k-Ni-J 

k-N+2 


h1 


k-N+iy 


k-N  +  1 


N-l„k-l 


y    -  nk_,  > 


l  N-  1     ,rk-N  +  l 

"k-N+iy 


k-l 


(4.37 


The  corresponding  coordinate  transformation  is  given  by 


dyA    (k-  N.N) 
civA 
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1 

0 

0 

0 

0 

-hk_N  +  ] 

1 

0 

0 

0 

~hk_N+1 

L    2 

_nk-N+2 

1 

0 

0 

N-2 


L  IN  — 2 

_nk-N  +  l 

i    N-l 
~nk-N  +  l 


.    N-2  i    N-2 

—  nk-N+2        _tlk-N+J 


2 

;-N+2 


.hN-l 


i    N-l 
fik-N+S 

~~ 0 


-h 


N-l 
k-1 


(4.38) 


These  definitions  are  sufficient  to  state  and  prove  the  theorems  presented 
in  the  remainder  of  this  chapter. 

2.     Generalized  Levinson  Algorithm 

The  Levinson  algorithm  (equation  (4.14))  can  be  extended  to  recursively 
compute  the  forward  and  backward  prediction  coefficients  defined  in  equations 
(4.26)  and  (4.31).  In  this  section  two  forms  of  the  algorithm  are  presented  and 
proved.  First  a  non-normalized  version  is  introduced,  then  using  this  result  a 
normalized  algorithm  is  developed. 

a.     Theorem  4.2:  Generalized  Levinson  Recursion  (Regular  Form) 

The  forward  and  backward  prediction  coefficients  defined  in 
equations  (4.26)  and  (4.31)  can  be  updated  using  the  following  recursive 
equations. 


(4.39a 


ihAN+ 

(k-1); 

=  !hAN(k+l)] 

-      ?Nrl 

M(k-N)j- 

1   -N     II 
<"k-l 

N 

~k-N  i     ! 

rN        1      1 

7  N+l 

(k-N)] 

=    hAN(k-N): 

k 

.  h,N(k-l)- 

rk-N 

nA 

\ 

1    ek-rl    ;      1 

\\ 


here 


k  ri-N    -N 

Pn  +  i   -    ^1ek  +  lrk-N 


(4.40) 


b.     Proof 

The    forward    and    backward    prediction    errors    are    scalars   so   their 
representation  is  identical  in  all  frames  of  reference.    Therefore,  we  can  write 
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■&  =  hAN(k+l)yA  (4.41a) 

=  KAN'(k+l)yv(k+l,N)  (4.41b) 
and 

rkN-N  =  hAN(k-N)yA  (4.42a) 

=  K^-(k-N)yA"(k-N,N)  (4.42b) 


where 


KAN(k  +  l)  =  hAN(k+l) — rr^ (4.43a) 

"l        ;           y        'ayA(k+l,N)  v           ' 

=  '0   ■  •  ■   0    -KkN_N+1    •  •  •  -KkN    1   0    •  •  •    0]                                              (4.43b) 

hAN(k.l)  =  K^(k+l)9yVft1,N)  (4-44) 

dy 


and 


KAN,-(k-N)  =  hAN(k-N) — ^ (4.45a) 

;  y3yA    (k-N.N)  V  ' 


=  [0    •  •.•    0   1    -K£N+1    ■  •  •    -Kjf    0    •  ■  •    oj  (4.45b) 

h*(k-N)  =  kA^(k-N)ayV>;N'N)  (4.46) 

dy 

The  normal  equations  in  the  primed  and  unprimed  coordinate 
systems  can  be  solved  for  KAN-(k-l)  and  KAN-(k-\).  This  is  once  again  (see  Section 
A  of  this  chapter)  a  straight  forward  matter  because  the  correlation  matrices. 
[E{yA'(k-4-i,N)y"'(k+l,N)}]  and  [E{yA"(k  N.N)y"  (k  S.S)}}.  will  contain  only 
diagonal  elements.    The  solutions  are 

KAN(k^l)  = 
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E|y(k-l)yk-w+1(k+l,N»              E<y(k+l)y*(k-l,N)> 
0    •  •  •    0    — ^-5 r-^-    •  "  •    — ^-j <~-    1    0    •  •  •    0 


E<Wyk-*  +  1(k+l,N))' 


E  (yk(k~l.N)): 


KAN-(k-N)  = 


(4.47) 


Eiy(k-N)yk-N+1(k-N,N)f  E<y(k-N)yk(k-N,N) 


0-01 


E<(yk-N"1(k-N,N))' 


E^(yk(k-N,N)): 


o-o 


(4.48) 


Consider  the  (k-N)-th  term  of  the  (N  +  l)  order  version  of  (4.47 


K£V(k^l) 


Eiy(k-l)yk--^k+l,N) 


:|(yk-N(k-l.N))4 


[4.49a) 


Ey(k-l)rk% 


E<(rkNN)2 


;4.49b) 


^k^V^N 


/ 


MrJ! 


k-  N 


(4.49c) 


eft, 


N 

rk-  N 


k 


(4.49d) 
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Similarly, 

KkN+V(k-N)  =  |  ;r7;  ;  P$+i  (4.50) 

:  ek+i !    1 

For   the    N  =  l    model    the    forward    error   prediction   coefficients    are 
(from  equation  (4.44)) 

[W(k+1)]  =  KAMk+l)  9yV(k+1»N)  (4.51a) 

dy 

=  [0   •     •    0    -K,,1    1   0   •  •  •   0]  (4.51b) 


I    I   e  0      |    I 

|0    ■  •  ■    o         '        '  :       Plk    1    0   •  •  •    0]  (4.51c' 

1  rk  I    I 


Equation  (4.39a)  yields 


0 


[hjflk+i)]  =  [bjP(k+l)]  -  />tk[hA°(k)j  7~~77  (4.52a) 

!     '   rk   I     I 

I     I   ft.0      II 

=  [0   ■  •  •   0   1   0   •  •  •   0]  -  Pik— — n  JO   •  •  •   0    1    0    •  •  ■    0  (4.52b) 

I   I  *k  I   I 

=  |0    •••   0    -P,k!.l|ek°;1.:,1      1    0    •••    :  (4.52c) 

!      I    rk    I      I 

Therefore,  the  recursion  of  equation  (4.39a)  holds  for  N=l.  We  assume  it  is  valid 
for  N  and  verify  it  for  N-j-1.    From  equation  (4.44)  we  have 

■hA»~>(k-l)    =  Kri(k,l)&yV(k;1-N|  (4.53a) 
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-KkN  + 

N  +  l 


Sf(k+l) 


Kkn_^+1(k+l) 


KkN+1(k+l)  +  h^Kjl^k+l)  + 


hkNN+1KkN_V(k+l) 


+  h.fK^k+l) 


(4.53b) 


=   ih^(k-l):   -   K^k+lJiO  0    1    -hkw_N+i 


-kN    0 


(4.53c) 


h,N(k-l) 


rk     N 


^,hAN(k-N) 


This  completes  the  proof  of  equation  (4.39a).  Vsing  similar 
arguments  the  backwards  recursion  of  equation  (4.39b)  can  be  verified. 

c.      Theorem  4.3:  Generalized  Levinson  Recursion  (Normalized  Form) 

The  normalized  prediction  error  coefficients  defined  in  equations 
(4.30)  and  (4.34)  can  be  recursively  updated  using  the  following  recurence 
relations; 


k-fl) 

k    N] 


=  e(pNk+1; 


a'(k-l) 
bJ(k-N) 


(4.54) 


where 


®(pLi 


Vi 


I P  N  + 1 , 


k 


/'N-l 


4.55' 


and  px+i  is  given  by  equation  (4.40). 
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d.     Proof 

The  proof  follows  directly  from  the  non-  normalized  version  of  the 
generalized  Levinson  algorithm  (equations  (4.39))  and  the  definitions  of  the 
normalized  prediction  error  coefficients  (equations    (4.29)  and  (4.33)). 

Using  (4.29)  in  (4.39a)  we  obtain 


aAN+1(k+l) 


ekN+1  |    !  aAN(k^l)  -  ,nk+1^    '  ^  '    1     !    I  rkN_N  |    |  bAN(k-N)]  (4.56a) 


ek  +  l     I  I     I    rk-N 


N 
ek  +  l  I  N/.  *  k      ,    N 


„N  +  1 

ek-i 


aAN(k  +  l)  -  ^+1bAN(k-N)]  (4.56b) 


Similarly,  using  (4.33)  in  (4.39b)  we  obtain 

*»aN+,(*-N)  =        ,'  ^J   !    [bAN(k-N)  -  „Nk+1aAN(k+l)]  (4.57 

I    !   rk-N   I     i 

The  proof  of  the  fact  that 


N  i     i       N       ,     i 

;k  +  l  !  !     I   rk-N  I  1 


e&i1         I  i  ^NV :  :      y/r-  (p^)2 


(4.58) 


is  relegated  to  the  next  section. 

The  initial  values  for  forward  and  backward  normalized 
autoregressive  parameters  are  obtained  by  setting  N=0  in  equations  (4.29)  and 
(4.33).    This  yields 

aA°(k+l)  =  JO   •  •  •    0        ,    |c1+] — -    0    •  •  ■    Oj  (4.59a) 

i  I  y 

bA°(k)  -    0   •  •  •    0   — -i 0   •  ■  ■    o  (4.59b) 

We  note  the  similarities  in  the  two  generalized  forms  of  the  Levinson 
algorithm  presented  in  this  section  with  that  presented  earlier  in  this  chapter.  A 
little  reflection  will  convince  that  equations  (4.14)  are  simply  a  special  case  of 
equations  (4.39). 
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3.     Error  Order  Update  Recursions 
a.     Theorem  4.4 

The  (N+1)  order  errors  can  be  calculated  from  the  N-th  order  errors 
through  the  recursion 


(4.60) 


where  0(p^+1)  is  defined  in  equation  (4.55). 
b.     Proof  (Outline) 

Using    the    normalized    form    of   the    Levinson    algorithm    (equation 
(4.54))  and  the  following  relationships 

ekN+1  =  af(k+l)/  (4.28b) 


-N  +  1 

-N 

«k  +  l 

Ck  +  1 

-N  +  1 
*-N 

=    ©(PN  +  l) 

-N 

rk-N 

rN  +  l   _       N  +  1, 
-k  +  1      -    aA         I 


k+l)y; 


rkNN  =  bAN(k-NV 


(4.32b) 


l£$  =  bAN+1(k-NV 

equation  (4.60)  is  easily  verified. 

We  now  return  to  the  proof  of  equation    (4.58).    That  is. 

I  ekN+,  !    |  '    I  -N  - 


rk-N 


i     i    „  N  + 1 
I    ek  +  l 


I    !  rkN  N 


Working  with  the  term  on  the  left-hand  side  of  the  equation 


1    I  eJ* , 

ck  +  l 


«k  +  l 


EUe? 


1     2 


E(ek^ 


.'  -  1  \  2  I 


(4.58) 


(4.61a) 


*>k+,ekN+i) 


T7T 


k  +  1  „  N  +  1 
k  +  1 


Eyk+1e 


(4.61b) 
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^N+iMy    rk-Nf ;     ek-i  i 


yk+1ekN+1 


T~1 


(4.61c) 


\/l  -  (/>n+i): 


(4.61d) 


Similar  arguments  can  be  used  to  verify  the  other  relationship  given  in  equation 
(4.58). 

Figure  4.3  symbolizes  a  single  stage  of  the  recursions  of  equation 
(4.60)  while  Figure  4.4  illustrates  a  third  order  analysis  model. 

The  interesting  feature  of  equations  (4.60)  is  that  they  do  not  make 
any  assumptions  about  the  nature  of  the  given  data.  The  data  values  need  not 
be  delayed  versions  of  each  other,  as  is  the  case  for  the  autoregressive  model. 
Any  set  of  data  values  can  be  used.  This  fact  will  be  significant  when  we  deal 
with  2-D  and  nonlinear  lattices. 

4.      The  Generalized  Schur  Algorithm 

In  this  section  an  algorithm  will  be  presented  which  allows  the 
calculation  of  the  partial  correlation  coefficients  in  a  direct  maimer.  It  will  b< 
shown  that  knowledge  of  the  correlation  matrix  is  sufficient  to  calculate  all  the 
reflection  factors  and  thus  solve  the  normal  equations  (by  use  of  the  Levinson 
algorithm).    Tin1  method  used  has  come  to  be  known  as  the  Schur  algorithm  [Ref. 

3]- 

In  order  to  obtain  the  desired  result  we  must  introduce  two  new- 
quantities,  defined  as 
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Figure  4.3:  Generalized  Lattice  Filter  Sections 
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Figure  4.4:  3-rd  Order  Generalized  Lattice  Filter 
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-N    „A 


aN^(k+l)  =  E^,y^=  a*(k-l)R 


./iA 


(4.62) 


/9NA+1(k-N)  =  E<kN_Ny4  =  b„N(k-N)R 


ma 


(4.63) 


where 


R"A  =  EJy»y 


(4.64) 


is  a  correlation  matrix. 

a.     Theorem  4.5:  Schur  Recursions 

The  quantities  a^(k-l)  and  /^(k-N)  defined  in  equations  (4.62)   and 
(4.63)  can  be  updated  according  to 


oNA+,(k+l) 
4i+i(k-N) 


9{pLi) 


aNA(k+l) 
^N(k-N) 


(4.65) 


and  the  partial  correlation  coefficient.  p^+u  can  be  calculated  from 
aNkN(k+l) 


Pn-i  = 


ak-N 


0]T*(k-N) 


[4.66) 


b.     Proof 

The  proof  of  equations  (4.65)  follows  directly  from  equations  (4.60) 
and  the  above  definitions  for  a^(k-  l)  and  $NA(k-N). 

Beginning  with  the  definition  of  the  partial  correlation  coefficient 
given  in  (4.40).  the  relationship  given  in  (4.66)  is  verified  as  follows 


L^ 
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E|ekN;,y4bAN(k-N) 


4.67a' 


[4.67b 


N  ,.i<-nLuN 


-    E^eir+Iyfc-N^_N(k-N) 


4.67c) 
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=  aNk-N(k+l)bk%(k-N)  (4.67d) 

But 

bk%(k-N)  =    :   :  yk1.N  |   |     =    [R(k_w)(1h_N)j1/2  (4.68) 

and 

/?Nk-N(k-N)  =  [R(x-N)(k-N)]i/2  (4_69) 

Therefore, 

/?Nk-*(k-N)  V          ' 

Using    the    initial    conditions    given  by    (4.59)    the    following    initial 
values  are  determined  for  the  Schur  algorithm 

R<k+1'A  (4.71a) 

[R(k+i)o     R(k  +  ,)i    .  .  .    R(k+i)K  (4.71b) 

(4.72a) 
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I  ;  1  i 

I  '9  A  1  k  V  - 

1     P 
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k  i    K 
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kA 


Rk0   Rkl    ■  ■       RkK  (4.72b) 

where  the  parenthesis  used  to  surround  the  first  indices  in  equation  (4.71)  simply 
indicate  that  they  are  fixed  at  the  value  indicated. 

The   Schur  algorithm  implies   a   filter  structure  identical   to   that    of 
Figure  4.4.    In  this  case  the  input  vectors  are  the  rows  of  the  correlation  matrix 
(normalized  by  the  square  root  of  the  diagonal  elements). 
5.      Synthesis  Model 

The  original  data.  Y.  can  be  synthesised  from  the  model  parameters 
obtained  from  the  Levinson  algorithm,  equation  (4.54).  It  is  also  possible  to 
regenerate   the   yA    directly    from   the    lattice    parameters.     The   desired   result    is 
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obtained  by  solving  equations  (4.60)  for  e^i,  and  rji"^- 

ek^l   =   V1   -    (Pn+i)     ek-l     +  *>N-lrk-N  (4.<oaj 


-N+l 


rk-N   =  n/1-  (/»Nk+i)V-N  -  PNk+1ekN+V  (4.73b) 

Equations  (4.73)  constitute  the  synthesis  model.  They  imply  a  structure 
similar  to  the  analysis  model  of  Figure  4.4.  but  with  the  direction  of  flow  of  the 
forward  error  signals  reversed.  The  processing  performed  at  each  stage  can  be 
visualized  as  in  Figure  4.5.  A  complete  third  order  synthesis  model  is  shown  in 
Figure  4.6.  Each  horizontal  path  in  Figure  4.6  represents  a  separate  synthesis 
model,  synthesising  a  different  component  of  Y.  The  coordinate  system  for  each 
of  these  models  depends  only  on  values  of  yA  which  appear  farther  down,  that  is 
they  have  a  smaller  value  of  the  index,  A. 

Compare  Figure  4.3b  and  4.5b.  It  is  apparent  that  the  behaviour  of  the 
backwards  error  signals  is  identical  in  the  two  cases.  Hence,  it  is  possible  to 
construct  a  synthesis  model  that  only  reverses  the  direction  of  the  forward  error 
corresponding  to  the  point  being  predicted.  This  assumes  knowledge  of  the  other 
inputs  (zero  order  forward  errors)  to  the  lattice.  Such  a  configuration  is 
illustrated  in  Figure  4.7. 

The  amount  of  knowledge  possessed  about  the  signals  used  for  the 
predictions  dictates  which  form  of  the  synthesis  model  should  be  used.  If  little  is 
known,  then  estimates  must  first  be  generated  which  can  then  be  used  in  the 
prediction.  This  corresponds  to  the  model  of  Figure  4.6.  If  complete  knowledge 
is  available  (either  from  initial  conditions  or  previous  predictions)  then  Figure  4.7 
can  be  used.  It  is  also  possible  to  construct  models  which  exploit  partial 
knowledge  of  the  input  signals  and  thus  fall  between  these  two  extremes.  In  this 
ease  the  known  signals  should  be  input  as  zero  order  forward  errors  while  the 
unknown  ones  must  be  estimated. 
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6.      Stochastic  Fourier  Series  Interpretation 

From  equations  (4.72)  and  (4.40)  we  can  deduce  the  following  result 

ek°+1  =  e^1  +   £  eLW-nW-a  (4.74a) 


N 
I 


e&t1*  E   I   I  ^+il   \pm*£-x  (4.74b) 


where  ek°+1  is  equivalent  to  yk+1.  These  expressions  offer  an  alternate  interpretation 
of  the  lattice  filter.  Equations  (4.74)  describe  a  stochastic  Fourier  series 
expansion  of  the  forward  error  sequence  where  the  basis  functions  are  given  by 
the  backwards  error  signals.  The  Fourier  coefficients  are  related  to  the  partial 
correlation  coefficients. 

This  concludes  the  review  of  existing  lattice  formulations.  In  the  next 
chapter  new.  multidimensional  extensions  of  this  theory  are  presented.  In 
Chapter  VI  these  results  are  used  to  derive  original  nonlinear  lattice  structures. 
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V.    TWO  DIMENSIONAL  LATTICE  STRUCTURES 

In  this  chapter  a  new  two-dimensional  lattice  structure  will  be  derived  and 
discussed.  Lattice  modelling  of  two  dimensional  fields  has  recently  received 
considerable  attention  [Ref.  43,44,45].  In  one  dimensional  lattice  modelling, 
updating  the  order  introduces  only  one  new  point  into  the  model  support.  An 
order  update  in  a  2-D  lattice  model  must  introduce  O(N)  new  points,  where  N  is 
the  model  order.  Several  different  solutions  have  been  suggested  to  this  problem. 
The  first  due  to  Marzetta  [Ref.  43,44],  uses  a  particular  ordering  of  the  data  to 
reduce  the  problem  to  one  dimension.  He  proposed  a  half-plane  support  which  is 
infinite  in  one  of  the  two  dimensions.  This  approach,  while  maintaining  several 
of  the  nice  characteristics  of  1-D  lattices,  such  as  correlation  matching  arid 
producing  a  minimum  phase  filter,  leads  to  very  long  delay  filters. 

A  different  approach,  proposed  by  Parker  and  Kayran  [Ref.  45]. 
simultaneously  introduces  many  points  into  the  support  when  the  model  order  is 
increased.  Their  filter  uses  a  quarter  plane  support  and  introduces  three 
parameters  at  each  order  update.  Therefore,  it  lacks  sufficient  parameters  to 
represent  all  classes  of  2-D  autoregressive  quarter  plane  filters.  More 
importantly,  it  lacks  the  property  of  orthogonality  so  that  the  cascading  of  stages 
does  not  lead  to  an  optimum  filter  (better  filters  are  possible  using  an  equivalent 
number  of  parameters).  It's  simplicity  is  attractive  and  good  results  have  been 
reported  using  this  approach  [Ref.  46]. 

The  theory  presented  here  maintains  features  of  both  previous  approaches.  It 
utilizes  the  generalized  lattice  theory  presented  in  Chapter  4  to  decompose  th( 
global  O(N)  point  update  into  0(N2)  single  point  local  updates.  It  maintains  the 
important  property  of  orthogonality  so  that  the  solution  at  all  stage?  is  optimum. 
Although  only  the  quarter  plane  support  case  is  presented  here  in  detail,  the 
theory  can  be  used  for  any  shaped  support.  It  is  shown  that  the  Levinson  and 
Schur  algorithms  (see  Chapter  4)  can  be  used  to  solve  the  2-D  linear  prediction 
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problem.  In  its  most  general  form  the  lattice  contains  0(N4)  parameters  while 
there  are  only  0(N2)  points  in  the  filter  support.  Several  structures  are  presented 
which  take  advantage  of  shift-invariance  and  reduce  this  requirement  to  0(\3). 

A.    GENERAL  FORM  OF  2-D  LATTICE  FILTER 

The  theory  used  is  exactly  that  presented  in  the  previous  chapter.  The  2-D 
structure  results  from  a  careful  selection  of  input  data.  To  illustrate  the 
proposed  2-D  lattice  structure  we  will  consider  a  2-D  linear  prediction  problem 
which  utilizes  a  quarter  plane  support.    The  2-D  data  field  is  given  by 

Y=/yV*=  y(XuX2)\  (5.1) 

where  (Ai,A2)  G  2LK  =   LK  x  LK 
where  LK  is  an  index  set  given  by 


LK  =  {o,  ...,K>  (5.2) 

Points  will  be  used  from  this  data  field  in  a   particular,  convenient   order.     We 
define  an  ordered  index  set 


2t  k 


6LK  =  <UO.O),(l,0),(0.1),(2.0),(0,2);  ....  (N,0),(0,N) 


(1;1),(2.1).(1.2).  ...,  (\.1),(1.N)..... 

(K-2,K),(K-1,K-  1),(K,K-1),(K-1,K),(K,K)|  (5.3) 

Other  orderings  are  possible  and  equally  valid.  This  one  is  chosen  merely  to 
illustrate  the  concepts.  The  desired  ordering  of  2LK  .  to  obtain  qLk  .  can  be 
accomplished  by  the  following,  computationally  efficient  algorithm 


108 


k  =  0 

for  mO  =  0  to  K 

for  nO  =  0  to  2*m0 
if  (mod(n0.2)  =  0) 
then  begin 
i  =  mO 
j  =  n0/2 
end 

else  begin 
i  =  n0/2 
j  =  mO 
end 
o2L*  (k)  =  2LK  (ij) 
k  =  k  +  1 
next  nO 
next  mO 

In  this  algorithm  qLk  (k)  has  been  used  to  describe  the  k-th  element  of  the  index 
set  qLk.  while  2LK(iJ)  has  been  used  to  indicate  the  (ij)-th  element  of  2LK.  The 
(K+l)2  elements  of  2LK  have  been  ordered  into  a  one  dimensional  index  set.  qLk. 
The  elements  of  qLk  can  be  numbered,  consecutively,  from  0  to  (K+l)2-l.  The 
notation  (k.l)  -  q  will  be  used  to  mean:  the  element  of  the  index  set  corresponding 
to  the  q-th  element  prior  to  the  element  (k.l).  For  example.  (2,0)-3  would 
indicate  the  element  (1.0)  (see  equation  (5.3)).  Occasionally  this  notation  will  be 
abbreviated  to  simply,  kl-q. 

Define  the  (q-1)  order,  normalized,  forward  error  associated  with  the 
prediction  of  the  point  y(k.l)  from  the  previous  (in  the  sense  of  qLk)  (q-1)  points, 
as 

ekr'  =  aA^(k,l)yV'  (5.4) 

where  the  implied  summation  over  (A,.A2)  e  2LK.  can  be  carried  out  in  any  order, 
as  long  as  all  components  are  considered.  It  is  preferred  to  think  of  (A,. A.,) 
belonging  to  2LK  rather  than  qLk  as  this  maintains  the  two-dimensional  character 
of  the  problem. 

The  a^fk.l)  can  be  interpreted  as  the  components  of  a  second  order 
covariant  tensor.    These  components,  for  a  range  of  indices  (A,.A2)  •'  (k.l)    q  (in  the 
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sense  of  (5.3)),  or  for  indices  (Aj.A2)  ">  (k-l)-  are  equal  to  zero.    When  (A,.A2)  =  (k.l). 
the  component 


ar^k,!) 


I    I  ekT 


•i  I    I 


(5.5) 


where  e^-1  is  an  unnormalized  version  of  ekT\ 

A  normalized  backward  error  associated  with  the  prediction  of  y((k,l)-q)  from 
the  next  q-1  points  of  qLk,  is  defined  by 

f£i  =  b^\((kA)-q)yX^  (5.6) 

As  with  the  forward  error  prediction  coefficients,  the  bAqV((k.l)-q)  can  be 
interpreted  as  components  of  a  second  order  covariant  tensor.  The  components 
b;\qV((k-l)-q)  equal  zero  for  the  range  of  indices  (A1;A2)  $  ((k.l)  — q)  or  (A,.A2)  >  (k,l). 
For  the  case    when    the    index  (A!,A2)  =  ((k,l)-q)  the  component 

l 


bk?:i((k.l)-q)  = 


nUi 


(5.7) 


where  r^lj  is  an  unnormalized  version  of  t^Z\. 
1.     Normalized  2-D  Levinson  Algorithm 
a.      Theorem  5.1 

The  2-D  prediction  error  coefficients  can  be  updated  in  order  using 
the  following  recursions 


a>V:(U) 


K 


©(O 


a^fk.l) 
b^l((k,l)-q) 


(5.8) 


where 


©(/»< 


an< 


l  -  (P?) 


] 


[5.9] 


pqM  -  Efer^i] 


(5.10) 
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b.     Proof  (Outline) 

The  proof  follows  directly  if  the  index  (k.l)  is  replaced  by  a  single 
index  which  runs  from  0  to  (K-l)2-  1  and  thus  indexes  the  elements  of  qLk.  The 
equations  (5.8)  are  identical  in  form  to  equations  (4.54)  and  the  proof  presented 
there  can  be  applied.  This  approach  is  equivalent  to  reordering  the  data  field 
into  a  vector  in  the  order  specified  by  qLk. 

An  alternate  proof  is  possible  by  generalizing  the  methods  of  Chapter 
4.  Alternate  coordinate  systems  could  be  introduced  (similar  to  (4.35)  and 
(4.37)).  The  necessary  transformations  can  be  found  and  all  the  steps  of  the 
proof  of  Theorem  4.2  can  be  generalized  to  these  higher  order  objects.  These 
concepts  will  not  be  explored  further  here  except  to  note  that  they  could  be 
extented  to  solve  problems  in  any  number  of  dimensions,  they  are  not  restricted 
to  the  two  dimensional  case  being  studied  here. 
2.  Normalized  2-D  Error  Order  Updates 
a.      Theorem  5.2 

The  two-dimensional  prediction  errors  can  be  updated  in  order 
according  to  the  following  equation 


Wl 


rkl-q 


ekT 

1 

hT- 

1 
q 

6Kkl)      _q_,  (5-11) 


b.     Proof  (Outline) 

The  proof  follows  from  equations  (5.8).    If  an  inner  product  is  formed 
on    both    sides   of  the    equation    (5.8a)    with    the    data   field    Y    equation    (5.11a) 
results.    Similar  arguments  can  be  employed  in  the  verification  of  (5.11b). 
3.      2-D  Form  of  Schur  Recursion 
Define  the  quantities 

vA 


aq_V(k.l)  =    E|ek1   V  2j  =   a.Vjk.lJR-'  "•  (5.12a) 

<W2((k,i)-q)  -  Etes:Jyvj  =  bAV4(w-q)R  W:A'  (5.12b) 


111 


w 


here 


A,A2AjA4  I    AjA2    A5A4 


:5.i3 


We  note  that  for  the  2-D  case  the  correlation  is  a  fourth  order  tensor. 

The  Schur  recursions  for  the  2-D  case  are  then  given  by  the  following 
theorem. 

a.     Theorem  5.3 

The  generalized  2-D  Schur  recursions  are  given  by 


*qV2(k.l) 

A,V2((M-q) 


=    ©(/>< 


"c-'hk,!) 


V.KMl-q) 


(5.14) 


and    p  kl  can  be  calculated  from 


A'q     = 


*qi 


k,I) 


5.15 


^((k.lj-q) 

b.     Proof  (Outline) 

The  proof  follows  identical  arguments  to  those  of  Theorem  4.5  and  so 
is  not  given  here. 

4.      2-D  Lattice  Structures 

The  derivations  presented  do  not  assume  shift-  invariance.  Models  are 
built  for  each  point  in  the  data  field  starting  with  the  point  y(K.K)  and  ending 
with  the  point  y(0.0).  All  the  models  are  not  equivalent,  however.  In  fact,  no 
two  models  are  identical.  The  only  model  that  is  quarter  plane  is  the  first  one. 
That  is  the  model  corresponding  to  the  point  y(K.K).  A  quarter  plane  model  for 
any  point,  y(m.n).  in  the  field  can  be  built  by  considering  an  appropriate  subset 
of  the  set  Y  (equation  (5.1)).  The  subset  would  sTarT  with  the  point  y(m-N.n-N) 
and  continue  in  a  quarter  plane  manner  until  the  point  y(m.n).  for  an  X-th 
(global)  order  model.    This  support  can  be  written 


D 


when 


y(A,.A2) 


:5.16) 
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(Ai,A2)  6    Lmn  —  Lm   x  Ln* 
where 

LmN  =  j(m-N),(m-N+l),  ....  ml  (5.17a) 

and 

LnN  =  |(n-N),(n-N+l),  ...,  n|  (5.17b) 

The  (N^l)2  terms  of  this  index  set  can  be  ordered  to  form 

02LmNn  =  |(m-N,n-N)I(m-N+lJn-N),(m-NJii-N+l))...1 

(m-lJn-l),(m,n-l),(m-l,n),(m,n)|  (5.18) 

The  subset  of  Y  given  by  (5.16)  and  the  ordering  of  (5.18)  are  illustrated  in 
Figure  5.1.  This  could  be  done  for  each  point  in  the  entire  data  field  Y.  yielding 
(K-l)2  models.  If  the  process  is  shift-invariant,  then  all  the  models  would  be 
identical.  Other  simplifications  in  the  model  are  possible  if  shift-invariance  can 
be  assumed.  These  will  be  the  topic  of  the  next  section.  In  addition,  if 
ergodicity  is  assumed  the  required  statistical  averages  can  be  replaced  by 
appropriate  spatial  ones. 

A  second  order  quarter  plane  2-D  lattice  model  for  generating  the 
prediction  error  associated  with  an  arbitrary  point  y(m.n)  is  illustrated  in  Figure 
5.2.  In  this  diagram  the  forward  and  backward  error  fields  are  indicated 
pictorially  rather  than  symbolically.  The  icons  used  are  defined  in  the  legend  on 
the  diagram.  The  large  squares,  at  each  stage,  indicate  the  entire  support  for  the 
second  order  model.  The  small  blank  squares  indicate  the  additional  support 
(besides  the  error  field  squares)  used  to  generate  the  given  error  fields.  For 
example,  the  forward  error  field  in  the  upper  right  hand  square  is  generated  by 
predicting  y(m.n)  from  all  the  remaining  data  points  in  the  large  square, 
including  the  one  indicated  as  a  backward  prediction  error. 
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The  doubly  hatched  squares,  corresponding  in  this  diagram  to  zero  order 
errors,  are  inputs  to  the  filter.  In  general,  (in  later  diagrams)  doubly  hatched 
squares  are  considered  to  be  inputs,  although,  they  may  not  always  correspond  to 
zero  order  errors. 

The  ordering  chosen  for  qLk  (equation  (5.3))  is  only  one  of  many  that 
could  have  been  chosen  to  implement  a  quarter  plane  model.  This  particular 
selection  was  made  for  ease  of  implementation  and  because  it  guaranteed  that  for 
every  2N+1  local  updates,  a  global  order  update  would  be  completed.  In  Figure 
5.2  the  filter  can  be  visualized  as  a  cascade  of  increasing  order  filter  sections.  For 
every  global  update.  0(N3)  local  order  updates  must  be  completed.  This  implies  a 
total  of  0(N4)  updates  for  an  N-th  order  filter.  In  general  this  is  too  large  a 
number  to  allow  these  filters  to  be  used  for  any  real  time  applications. 

In  the  next  section  the  problem  of  reducing  the  complexity  of  the  2-D 
lattice  filter  is  examined  and  some  solutions  are  proposed. 

B.    REDUCED  COMPLEXITY  2-D  LATTICE  FILTERS 

The  assumption  of  shift-invariance  allows  certain  of  the  backward  prediction 
errors  to  be  considered  to  be  shifted  versions  of  each  other.  This  eliminates  some 
calculations.  Various  structures  are  possible  depending  on  the  types  of  shifts 
introduced.  We  note  that  no  single  type  of  shift  (neither  z,-1,  z2_1,  z,_1z2-1)  will 
introduce  all  the  new  data  elements  into  the  support  that  are  required  for  a 
global  order  update.  Because  of  this,  additional  prediction  errors  will  have  to  be 
introduced  at  each  stage.  This  reduces  somewhat  the  advantage  gained  by  the 
shift -in  variance  assumption. 

Two  types  of  delays  will  be  considered  in  detail.  Initially,  several  models 
involving  diagonal  shifts  are  examined.  Later,  a  model  involving  a  horizontal 
shift  is  discussed.  We  begin  by  introducing  a  diagonal  shift  operator.  D.  This  is 
equivalent  To  multiplication  by  z,  'z2  '  in  the  bivariate  z-transform  domain. 

Because  of  the  assumed  shift  invariance  the  following  statements  can  be 
made 
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R(m.n.m-Ln-j)  =  E<y(m,n)y(m-i,n-j)|'  (5.19a) 

=  B&{y(m>n)]D[y(m-i,n-j)lJ  (5.19b) 

=  E|y(m-l,n-l)y(m-i-l,n-j-l)|  (5.19c) 

=  R(ij)  (5.19d) 

where  R()  is  a  correlation  function.  The  correlation  is  only  a  function  of  the 
relative  positions  of  the  two  points  not  their  absolute  positions.  If  we  adopt  for  a 
moment  the  Hilbert  space  formulation  (see  Appendix  A)  we  conclude  that  the 
diagonal  shift  operator  is  an  inner-product  preserving  operator  and  so  the  use  of 
the  shifted  versions  of  the  backward  error  signals  is  justified. 

Consider  the  structure  illustrated  in  Figure  5.3.  It  represents  a  third  order 
quarter  plane  lattice  filter.  At  each  global  stage,  (2N-l):-3  lattice  coefficients 
are  introduced.  Therefore,  an  N-th  order  model  requires  0(N3)  parameters. 
Notice  thaT  at  each  stage  two  new  errors  are  introduced.  They  each  require  the 
solution  of  an  (N-l)2  point  prediction  problem.  For  small  N  this  is  an 
insignificant  number,  however,  for  large  N  it  becomes  overwhelming  and  the 
required  number  of  parameters  again  becomes  0(N4).  It  is  difficult  to  analyze  this 
structure  analytically,  as  the  index  sets  for  each  prediction  error  are  different. 
The  support  for  different  errors  follow  different  patterns.  This,  and  the 
complexity  issue  make  it  a  structure  that  is  really  only  of  academic  interest. 

The  structure  of  Figure  5.4  is  a  True  O(v')  parameter  model.  It  avoids  The 
addiTion  of  The  new  error  signals  aT  each  sTage  by  introducing  them  aT  The  ouT^et. 
The  structure  has  a  sup)port  thaT  differs  slightly  from  quarter  plane.  A  more 
signifiranT  drawback,  however,  is  ThaT  the  maximum  order  of  The  filter  must  be 
fixed  at  the  start.  If  the  maximum  order  is  overestimated  then  some  unnecessary 
computations  will  have   been   performed.     If  on   the   other   hand,  the  maximum 
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Figure  5.3:  Reduced  Complexity,  2-D  Quarter  Plane  Lattice  Filter. 
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required  order  was  underestimated,  then  a  great  price  must  be  paid  to  increase  it. 
However,  this  is  considered  to  be  a  superior  structure  because  of  the  regularity 
and  complexity  reduction  it  offers.  Therefore,  a  more  detailed  analysis  of  this 
model  will  be  performed. 

The  ordered  index  set  (equivalent  to  (5.18))   in  this  case  (for  N=2)  is  given 
by 

oLm2n  =  |(m-2,n-2),(m-l,n-3),(m-3,n-l),(m-l,n-2),(m-2,n-l), 

(m-l.n-l).(m.n-2),(m-2.n).(m,n-l),(m-l,n),(m.n)l  (5.20) 

The  following  relations  can  then  be  used  to  advantage  to  update  the  backwards 
errors 

r(°m_x,„-i)=  D,r(0m,n);  (5.21a 

Fi-2,n-i)=  D[r£,-id  (5-21b 

H2m-l.n-2)    =     Dir^.n-!)]  (5.21c 

r(m-3.n-i)  =  D!r(sm_2.n)]  (5.21d 

Vm-l.n-S)    =     D^.n-2);  (5-21e 

-5  t\  -5 


r(m-'.'.n-2) 


Di,.-,,,  (5.2lf 


In  general,  whenever  (m-i.n-j)  equals  (m.n)-q.  then 


—  q 

r(m -  l-  l.n-  j-  1) 


(5.22 


Equation   (5.22)   is  a  simple  rule  for  exploiting  the  shift-   invariance  of  the  data 
field. 

One  final  reduced  complexity  lattice  model  will  be  introduced.  It  will  serve 
to  illustrate  the  variety  of  structures  possible  and  in  particular  will  yield  a  model 
for  which  an  especially  convenient  synthesis  model  can  be  constructed.  The 
model  incorporates  a  horizontal  (z2  ')  shift,  rather  than  the  diagonal  shift  used  in 
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Figure  5.4:  Reduced  Complexity.  0(NS),  2-D  Quarter  Plane  Lattice  Filter 
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the    previous    models.     A    horizontal    shift    also    is    a    inner    product    preserving 
mapping  so  that  its  use  is  permitted. 

The  ordered  index  set  used  for  this  model  is  given  by  (for  the  second  order 
case) 

02L^n  =  |(m-2,n-2),(m-l,n-2),(m,n-2),(m-2,n-l), 

(m-l,n-l)!(m,n-l),(m-2,n),(m-l,n),(m,n)i  (5.23) 

Define  a  horizontal  shift  operator,  H,  by  the  relation 

y(m.n-l)  -  H(y(ra,n>]  (5.24) 

The  following  relations  hold  due  the  assumed  shift-  invariance  of  the  given  data 
field 

F(°m,n-i)=  Hlr(°m,n)]  (5.25a) 

r(m-i,n-i)  =  Hr(ni^ln)  (5.25b) 

r(m-2.n-l)   =    HJF(m_2ln)  (5.25c) 

r"(m-s,n-i)  =  H[r(m_SiI1)  (5.25d) 

r(m-4,n-l)    =     H  ,r(m-4.n).  (5.25e) 

Hm-  S.n-1)   -    Hf|l-5,r,|  (5.25f) 

In  general,  whenever,  (m-i,n-j)  equals  (m,n)-q  then 

Hm-i.n    j-1)    :      Hr,^,,,,,.  j)  (5.26) 

Using  these  simplifying  relations  and  equations  (5.11).  the  model  of  Figure 
5.5  can  be  deduced.  This  model  still  contains  O(N')  parameters,  however,  the 
actual  number  is  only  about  one  quarter  of  that  required  by  the  previous 
structure  (Figure  5.4.) 


121 


This  algorithm  shares  with  the  previous  algorithm  the  shortcoming  of  having 
to  estimate  the  maximum  order  of  the  filter  at  the  outset.  Despite  this,  it  is 
believed  that  this  model  offers  a  good  compromise  (probably  the  best  to  date) 
between  model  accuracy  and  implementation  complexity.  In  the  next  section  it  is 
demonstrated  that  the  synthesis  form  of  this  algorithm  has  some  particularly 
desirable  properties.  In  Chapter  7,  it  will  be  shown  that  this  algorithm  is  well 
suited  for  highly  parallel  VLSI  implementation. 

C.    SYNTHESIS  MODEL 

The  synthesis  results  of  the  previous  chapter  are  easily  extended  to  the  2-D 
case  being  considered.  The  data  fields  can  be  regenerated  from  a  knowledge  of 
the  lattice  coefficients  and  the  forward  error  fields,  it  is  not  necessary  to  explicitly 
calculate  the  forward  and  backward  error  prediction  coefficients. 

The  desired  result  is  obtained  by  solving  equation  (5.11)  for  e^"1  and  rk1_q. 
This  yields 


ekr2  =  n/1-   (O'eiH  +  pfrZll  (5.27a) 


r^q  =  v/i-  (/»qu)2n3=S  -  p?*u  (5-27b) 

These  equations  establish  the  method  for  regenerating  the  original  data  field. 
They  describe  the  processing  that  must  be  carried  out  at  each  stage  of  the 
synthesis  process.  As  an  example  of  their  application,  consider  the  second  order 
synthesis  model  pictured  in  Figure  5.6.  It  is  the  synthesis  counterpart  of  the 
reduced  complexity  analysis  model  of  Figure  5.5.  In  order  to  regenerate  the 
original  data  field  processing  should  be  carried  out  horizontally  by  rows.  For  this 
second  order  model,  two  rows  and  two  columns  of  initial  conditions  must  be 
specified.  The  required  zero  order  forward  error  sequences  will  always  be 
available  from  either  the  given  initial  conditions  or  from  previous  estimates. 

It  will  be  shown  that  for  VLSI  implementations  it  will  be  more  convenient  to 
estimate  all  the  zero  order  forward  error  sequences.  This  necessitates  that  three 
residuals  to  be  input  and  that  all  forward  error  channels  be  reversed.  Such  a 
structure  is  illustrated  in  Figure  5.7. 
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A  synthesis  model  corresponding  to  the  analysis  model  of  Figure  5.4  would 
have  to  process  the  data  along  diagonals.  The  required  zero  order  forward  error 
sequences  would  not  be  available  at  each  stage  in  this  case  and  so  only  a  model 
analogous  to  the  one  of  Figure  5.7  can  be  used.  This  would  require  the  provision 
of  five  input  residual  signals. 

D.    SYSTOLIC  IMPLEMENTATIONS 

In  the  past,  most  signal  processing  algorithms  were  implemented  largely  in 
software  due  to  their  high  complexity.  More  recently,  with  the  advent  of  VLSI 
technology,  there  has  been  a  shift  towards  specialized  hardware  implementations. 
These  offer  increased  performance  at  a  low  per-unit  cost.  A  particularly 
promising  class  of  implementations,  first  suggested  by  Kung  and  Leiserson  [Ref. 
49].  are  the  so  called  systolic  arrays.  These  attempt  to  partition  the  required 
computations  in  time  and  space  over  an  array  of  identical  processing  elements,  in 
order  to  increase  throughput. 

Kung  [Ref.  50:  pp.  869]  defines  a  systolic  array  as  a  network  possessing  the 
following  features: 

a)Svnchrony:  The  data  are  rhythmically  computed  (timed  by  a  global  clock) 
and  passed  through  the  network. 

b)Regularity  (i.e..  Modularity  and  Local  Interconnections):  The  array  should 
consist  of  modular  processing  units  with  regular  and  (spatially)  local  intercon- 
nections.   Moreover,  the  computing  network  may  be  extended  indefinitely. 

c)Ten:poral  Locality:  There  will  be  at  least  one  unit-time  delay  allotted  so  that 
signal  transactions  from  one  node  to  the  next  can  be  completed. 

d)  Pipelinabilit\  (i.e..  O(M)  Execution-Time  Speed-l  p):  A  good  measure  for  the 
efficiency  of  the  array  is  the  following 


Speed-  Up  Factor 


Processing  Time  in  a  Single  Processor 
Processing  Time  m  the  Arra\   Processor 


A  systolic  arrav  should  exhibit^  a  linear-rate  pipelinability.  i.e..  it  should  achieve 
an  O(M)  speed-up.  in  term'-  of  processing  fates,  when  \"1  is  th<  number  of  pro- 
cessor elements  (PE"s). 


Methods  have  been  proposed  for  transforming  algorithms  into  Systolic 
implementations  beginning  with  either  an  algorithmic  description  (see 
Moldovan  [Ref.  51])  or  a  signal-flow-graph  (SFG)  description  (see  Kung  [Ref. 
50]).    In  this  chapter  we  shall  be  using  the  second  of  these  methods  to  transform 
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one  of  the  2-D  lattice  structures  into  systolic  form.  Rather  than  present  a 
detailed  discussion  of  the  rules  used  in  this  method,  we  shall  simply  state  them 
and  then  apply  them  to  produce  the  desired  systolic  array.  It  is  hoped  that  this 
example  will  clarify  the  procedures  used.  If  further  insight  is  desired  the  reader  is 
referred  to  the  cited  reference. 

1.  SFG  Transformation  Procedure 

The  procedure  used  is  based  on  a  cut-set  approach.  According  to  Kung 
[Ref.  50,  pp.  870]  a  cut-set  is  defined  as: 

A  cut-set  in  an  SFG  is  a  minimal  set  of  edges  which  partitions  the  SFG  into  two 
parts. 

He  proposes  and  proves  that  the  following  two  rules  that  can  be  used  to 
transform  any  computable  SFG  into  a  systolic  array: 

Rule  (i)  Time-Scaling:  All  delays  D  may  be  scaled,  i.e.,  D  —  D'.  by  a  single  po- 
sitive integer  .  Correspondingly,  the  input  and  output  rates  also  have  to  be 
scaled  by  a  factor     (with  respect  to  the  new  time  unit  D')  ...  . 

Rule  (ii)  Delay-Transfer:  Given  any  cut-set  of  the  SFG,  we  can  group  the  edges 
of  the  cut-set_  into  "inbound  edges'f  and_ "outbound  edges",  depending  upon  the 
directions  assigned  to  the  edges.  Rule  (ii)  allows  advancing  k  (D')  time  units  on 
all  the  outbound  edges  and  delaying  k  time  units  on  the  inbound  edges,  and  vice 
versa.  It  is  clear  that,  for  a  (time-invariant)  SFG.  the  general  system  behaviour 
is  not  affected  because  the  effects  of  the  lags  and  advances  cancel  each  other  in 
the  overall  timing.  Note  that  the  input-input  and  input-output  timing  relation- 
ship? will  also  remain  exactly  the  same  only  if  they  are  located  on  the  same  side. 
Otherwise  they  should  be  adjusted  by  a  lag  of  — ktime  units  or  an  advance  of -k 
time  units. 

2.  Systolic  Implementation  of  2-D  Lattice  Filter 

Using  the  two  rules  given  in  the  previous  section  the  2-D  lattice  synthesis 
model  of  Figure  5.7  will  be  mapped  into  a  systolic  array.  There  is  some 
flexibility  in  the  design,  the  result  not  being  unique.  The  first  choice  that  must 
be  made  is  that  of  the  operation  that  is  to  be  performed  by  the  basic  PE.  In  the 
case  being  considered  several  convenient  choices  are  possible.  The  simplest 
element  would  be  a  multiplier-adder.  A  slightly  higher  level  operation  would  be 
that  of  the  >ingle  lattice  >ection  given  by  Figure  4.5.  A  still  higher  level  operation 
is  conceivable  by  grouping  several  of  the  lattice  sections  together.  We  -hall  use 
the  second  of  these  choices  as  it  illustrates  the  general  procedure  without  the 
added  complexity  inherent  in  the  lower  level  implementation. 
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The  mapping  will  be  done  in  stages.  Initially,  the  diagram  of  Figure  5.7 
will  be  redrawn  in  an  SFG  format.  Rule  (i)  will  be  used  to  scale  all  the  delays 
appropriately.  Then,  by  successive  application  of  Rule  (ii).  the  SFG  will  be 
temporally  localized  (it  is  already  spatially  local.)  The  steps  used  are  outlined  as 
follows; 

(1)  In  Figure  5.8  the  algorithm  of  Figure  5.7  has  been  redrawn  in  a  SFG  form. 
In  this  and  subsequent  diagrams  delays  are  indicated  by  the  letter  D.  The 
lattice  section  of  Figure  4.5,  as  before,  is  indicated  by  the  circles  at  the 
nodes. 

(2)  Using  rule(i),  all  delays  are  scaled  by  a  factor  of  6.  That  is,  D  —  6D'.  This  is 
indicated  in  Figure  5.9.  The  input  and  output  signal  rates  must  also  be 
scaled  by  this  factor  of  6. 

(3)  Using  the  cut-sets  indicated  in  Figure  5.10  the  delays  are  redistributed  so 
that  temporal  locality  is  achieved.  The  resulting  SFG  is  indicated  in  Figure 
5.11.  Until  now  the  processing  at  each  node  was  assumed  not  to  take  any 
time.  The  delays  going  into  each  node  can  be  combined  with  the  lattice 
sections  and  be  used  to  account  for  the  processing  time.  In  this  way.  the 
structure  will  appear  as  in  Figure  5.12.  In  this  last  figure,  the  nodes  have 
been  shaded  to  indicate  that  the  operation  being  performed  within  them 
consumes  one  time  unit. 

3.      Additional  Remarks 

In  this  section  we  have  shown  that  the  2-D  Lattice  structures  derived  in 
this  chapter  are  amenable  to  a  systolic  implementation.  This  is  significant  as  the 
processing  of  2-D  data  fields  such  as  images  in  real-time  requires  high  data  rates. 
These  rates  can  only  be  achieved  in  practice  through  the  use  of  super-computers 
or  specialized  hardware.  Due  to  the  high  cost  of  super-computers  the  second 
alternative  is  the  more  practical.  With  the  costs  of  VLSI  production  rapidly 
decreasing,  it  is  now  cost  effective  to  produce  dedicated  chips  even  in  very  small 
quantities.  For  large  scale  productions  the  cost  can  be  amortized  over  a  large 
number  of  chips,  yielding  a  low  per-unit  cost. 
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Although  only  a  single  specific  implementation  has  been  presented  here, 
it  indicates  the  ease  with  which  the  other  algorithms  discussed  in  this  thesis  may 
be  transformed  into  forms  which  can  be  efficiently  implemented  in  silicon. 

E.    SIMULATION  RESULTS 

The  theory  has  been  proven  by  computer  simulation.  Two  different  order 
models  were  excited  by  a  white  noise  process.  Several  different  order  estimates  of 
each  spectrum  were  generated  and  compared  to  the  originals. 

1.  Example  5.1 

The  first  model  simulated  was  described  by 

y(m,n)  =   .295y(m-l,n)  -  .470y(m,n-l)  +  0.0y(m-l,n- 1) 

-  .055y(rn-2,n)  +  .007y(m,n-2)  +  .003y(m-2,n- 1) 

-  .015y(m-l,n-2)  +  .022y(m-2,n-2)  -  u(m,n)  (5.28) 

where  u(m,n)  was  a  2-D  .zero-mean  white  noise  process. 

The  original  spectrum  is  illustrated  in  Figure  5.13  while  the  first,  second, 
third,  and  fourth  order  estimates  are  given  in  Figure  5.14.  Notice  that  the 
original  model  is  only  second  order  so  that  the  third  and  fourth  order  estimates 
can  be  used  to  examine  the  effects  of  over  modelling.  As  can  be  seen  from  the 
figures  the  estimated  spectra  correspond  very  closely  to  the  originals.  Over 
modelling  did  not  noticeably  degrade  the  accuracy  of  the  estimates. 

The  actual  algorithm  used  was  that  of  Figure  5.2.  Although  it  is 
unnecessarily  complex,  it  is  the  most  straight  forward  to  implement.  The 
generalized  2-D  Schur  algorithm  was  used  to  generate  all  the  required  lattice 
parameters.  The  computer  subroutines  used  to  accomplish  this  simulation  are 
included  in  Appendix  B. 

2.  Example  5.2 

The  second  simulation  used  the  following  higher  order  a.utoregressive 
equation  to  generate  the  data 

y(m,n)  =   -.47y(m-l.n)  -  .03y(m,n-  1)  +  .195y(m  -  l.n  -1) 
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Figure  5.8:  2-D  Lattice  Filter  of  Figure  5.7  in  SFG  Form. 
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Figure  5.9:  The  Delays  are  Scaled  by  a  Factor  of  Six. 


131 


CU 


IS.  U> 


u> 


00 

to 
in 


«■-'*- 


-r1 


o — "^— <» 


CO 


in 


'L-iJ'. 


o  ^  o  «q  o 


t.-j,---^-^ 


co 


,_„._ 


(n _...,,.._, -,,.- 


_  X 


H- 


f -T 


J_'L 


4 


T 


— » — y 


—  i — t 

i 

± 


CU 


00 
CO 

in 


CO 

in 


CO 

in 


Figure  5.10:  SFG  for  2-D  Lattice  Filter  Showing  Cut-Sets  to  be  Used 
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Figure  5.13:  Original  Spectrum  For  Example  5.1. 
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a.    First  Order  Lattice  Approximation  of  Spectrum  of  Example  5.1. 


b.    Second  Order  Lattice  Approximation  of  Spectrum  of  Example  5.1. 
Figure  5.14:  Lattice  Approximations  Of  Spectrum  of  Example  5.1 
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c.    Third  Order  Lattice  Approximation  of  Spectrum  of  Example  5.1. 


d.    Fourth  Order  Lattice  Approximation  of  Spectrum  of  Example  5.1. 
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-  .015y(m,n-2)  +  .055y(m-l,n-2)  -   .003y(m-2,n-2) 
+  .0067y(m-3,n)  -  .015y(m.n-3)  -  .022y(m-2.n-  3) 
-+  .033y(m,n-4)  -  .085y(m-l,n-4)  -  .002y(m-4,n-2) 

-  .0001y(m-2,n-4)  +  .0001y(m-4,n-4)  +  u(m,n)  (5.29) 

where  u(m,n)  is  a  zero-mean  white  noise  process. 

The  spectrum  of  this  model  is  shown  in  Figure  5.15.  The  first  through 
fourth  order  estimates  of  the  spectrum  are  shown  in  Figures  5.16.  The  general 
shape  of  the  spectrum  is  identified  in  the  first  order  model,  although  the  fine 
detail  is  not  introduced  until  the  fourth  order  model.  The  position  and  relative 
magnitude  of  the  peaks  in  the  spectra  are  identified  with  great  accuracy  in  the 
fourth  order  estimate. 

In  the  next  chapter  lattice  models  are  applied  to  the  solution  of  the 
autoregressive  nonlinear  modelling  problem. 
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Figure  5.15:  Original  Spectrum  For  Example  5.2. 
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a.    First  Order  Lattice  Approximation  of  Spectrum  01  Example  5.2. 


b.    Second  Order  Lattice  Approximation  of  Spectrum  of  Example  5.2. 
Figure  5.16:  Lattice  Approximations  Of  Spectrum  of  Example  5.2 
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c.    Third  Order  Lattice  Approximation  of  Spectrum  of  Example  5.2. 


d.    Fourth  Order  Lattice  Approximation  of  Spectrum  of  Example  5.2. 
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VI.    NONLINEAR  LATTICE  STRUCTURES 

In  the  previous  two  chapters  lattice  structures  have  been  examined  in  great 
length  and  have  been  applied  to  the  problem  of  2-D  autoregressive  modelling.  In 
this  chapter  we  return  to  the  problem  of  nonlinear  autoregressive  modelling  with 
the  hope  of  solving  the  problem  using  a  lattice  structure. 

There  have  been  several  attempts  in  the  literature  to  fit  lattice  filters  to  a 
Volterra  like  model.  The  first  was  due  to  Parker  and  Perry  [Ref.  15].  They 
proposed  a  multichannel  lattice  implementation  of  an  autoregressive-moving 
average  nonlinear  model.  Their  model  was  capable  of  providing  an  update  in  time 
of  some  terms  of  the  expansion.  It.  however,  does  not  introduce  all  the  terms  that 
are  needed  to  constitute  the  next  higher  order  model. 

A  recent  proposal  by  Zarzycki  and  Dewilde  [Ref.  17]  is  a  true  generalization 
of  the  Levinson  algorithm  to  the  case  of  the  autoregressive  Volterra  type  model. 
In  their  work  they  considered  a  non-stationary  nonlinear  structure  and  showed 
that  the  stationary  and  linear  models  can  be  treated  as  special  cases.  Their 
model  provides  a  true  update  in  time  order,  introducing  all  the  required  terms. 
They  found,  as  we  did  with  the  2-D  model,  that  not  all  terms  could  be 
recursively  generated  and  that  some  would  have  to  be  introduced  at  each  stage 
by  other  means.  For  these  non-linear  models  only  one  error  signal  must  be 
injected  at  each  stage  not  two  as  in  the  2-D  case  (if  a  triangular  kernel  is  used.) 

The  model  proposed  in  this  chapter  is  based  on  the  alternate  tensor  form  of 
the  nonlinear  model  developed  in  Chapter  3.  Recall  that  the  model  was  based  on 
interchanging  the  roles  played  by  time  order  and  nonlinear  order.  The  lattice 
model  presented  here  thus  becomes  recursive  in  nonlinear  order.  This  means,  for 
example,  that  the  optimum  cubic  model  is  calculated  from  a  knowledge  of  the 
optimum  quadratic  model.  The  theory  is  based  on  the  generalized  lattice 
concepts  of  Chapter  4  and  the  2-D  lattice  structures  developed  in  Chapter  5. 
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Although  only  the  nonlinear  order  update  model  is  discussed,  some  reflection  will 
show  that  a  time  update  can  also  be  performed  using  the  proposed  model. 

For  simplicity  in  illustrating  the  concepts  only  models  involving  two  delays 
will  be  considered  in  this  chapter.  The  theory  is  equally  valid  and  easily  extended 
to  more  complex  situations. 

A.    GENERAL  NONLINEAR  LATTICE  MODEL 

As  with  the  2-D  model  of  the  previous  chapter,  the  theory  used  is  that  of  the 
generalized  error  order  updates  given  in  Chapter  4.  A  careful  choice  of  the  input 
data  will  yield  the  desired  results.  We  begin  by  summarizing  briefly  the 
nonlinear  model  to  be  used,  for  the  case  of  two  delays. 

The  estimate  of  the  model  output  is  given  by  (3.53) 


Av, 


(6.i; 


y(k)  =  y   "(k     1)    •■  ■    y  N(k-N)HAi    . .  ,N 
To  avoid  unnecessarily  complex  algebra  we  consider  the  case  when  N=2, 

y(k)  =  yAl(k-l)yAj(k-2)HV2  (6.2) 

The    tensor    product    y  '(k-  l)y  2(k-2)    defines    a   second   order   tensor,    Y(k),   with 
components  given  by 

Y(k)=  iyAl(k-l)yA2(k-2)]  = 


1  y(k-2)  yM(k-2) 

y(k-l)        y(k-l)y(k-2)        y(k-  l)y(2)(k-  2) 

y,2)(k-l)    y(2)(k-  l)y(k-2)    y,2)(k- 1  )y(2)(k- 2) 


ylp;(k     ])    y(Fl(k     l)y(k-  2)    y^(k-l)y(2>(k     2] 


yW(k-2) 
y(k-l)yW(k-2) 

y(2)(k-l)y(p)(k-2) 


(p,(k     l)v(pl(k     2) 


(C.3) 

This  tensor  can  be  considered  to  be  a  shift-varying  data  field  and  can  thus  be 
modelled  using  the  2-D  lattice  model  of  Figure  5.2.  A  tensor  of  the  form  given  in 
(6.3)  can  be  formed  for  each  time.  k.    If  the  process  is  time-varying  then  each  of 
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these  tensors  must  be  separately  modelled.  If  time-invariance  can  be  assumed, 
then  only  one  model  need  be  developed.  In  that  case  the  required  correlations 
can  be  calculated  over  the  ensemble  of  such  tensors. 

A  shift  along  the  diagonal  (or  horizontal  or  vertical  )  of  the  Y(k)  tensor  is  not 
a  inner-product  preserving  operation.  This  prevents  the  application  of  the 
reduced  complexity  algorithms  of  Chapter  5  to  this  nonlinear  problem. 

The  2-D  lattice  model  of  the  tensor  Y(k)  does  not  offer  a  complete  solution. 
Notice  that  y(k),  the  sample  that  is  to  be  estimated,  does  not  appear  in  equation 
(6.3).  Two  possible  solutions  to  this  problem  exist.  First,  the  ordered  index  set. 
(5.3),  can  be  extended  by  one  element  so  that  y(k)  is  included  in  the  model.  This 
has  the  effect  of  adding  a  channel  to  the  general  2-D  lattice  structure  of  Figure 
5.2.  A  conceptually  different  solution  is  to  first  model  the  tensor  Y(k)  using  the 
results  of  the  previous  chapter.  The  backwards  error  signals  that  result  can  be 
used  as  a  basis  for  a  Fourier  expansion  of  the  the  sample  y(k)  (see  equation 
(4.73)  for  details.)  Both  of  these  approaches  lead  to  identical  models  but  do 
provide  alternate  interpretations.  The  second  approach  will  be  the  one  used  in 
the  derivations  of  this  chapter  as  it  allows  the  results  of  the  previous  chapter  to 
be  applied  with  little  modification. 

The  tensor.  Y(k).  can  be  considered  to  be  a  two-  dimensional  data  field  given 
by 


Y(k)  =  y »(k-l)y  :(k-2)j  (6.4] 

whore  (A,.A2)       2V        lp    ■   U 
where  Lp  is  an  index  set  given  by 

....  pi  (6.5) 


Points  will  be   used  from   this  data   field   in   a   particular,  convenient   order.     We 
define  an  ordered  index  set 
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&p  =  |(p,p),(p-l,p),(p,p-l),  ...,  (0,p).(p.O), 

(p-l,p-l),(p-2,p-l),(p-l,p-2).  ...,  (0.p-l);(p-1.0), 

...,(1,1),  (0,1),(1,0),  (0,0)1  (6.6) 

The  (p-t-l)2  elements  of  2LP  have  been  ordered  into  a  one  dimensional  index 
set,  qLp.  The  elements    of  oLp  can  be  numbered,  consecutively,  from  0  to  (p^l)2-l. 

As  in  Chapter  5.  the  notation  (m,n)  -  q  will  be  used  to  mean;  the  element  of 
the  index  set  corresponding  to  the  q-th  element  prior  to  the  element  (m,n).  For 
example.  (1.0)  -  2  would  indicate  the  element  (1,1)  (see  equation  (6.6)).  This 
notation  will  often  be  abbreviated  to  simply  mn-q. 

Define  the  (q-1)  order,  normalized,  forward  error  associated  with  the 
prediction  of  the  element  ym(k-l)yn(k-2)  from  the  previous  (in  the  sense  of  qLp) 
(q-1)  points,  as 


e^1 


ajftfm.njy  '(k-l)y  2(k-2)  (6.7) 


where  the  implied  summation  over-  (Ai,A2)  ~  2LP,  can  be  carried  out  in  any  order, 
as  long  as  all  components  are  considered.  It  is  preferred  to  think  of  (Ai,A2) 
belonging  to  2LP  rather  than  £LP  as  this  maintains  the  multi-dimensional 
character  of  the  problem. 

The  ay^Jjm.n)  can  be  interpreted  as  the  components  of  a  second  order 
covariant  tensor.  These  components,  for  a  range  of  indices  (Aj,A2)  '  mn-q  (in  the 
sense  of  (6.6)).  or  for  indices  (A,,A2)  >  (m,n)  are  equal  to  zero.  For  the  case  when 
(A,.Ao)  =   (m.nj.  the  component 


amVfm.n; 


(6.s; 


where  e^,,1  is  an  unnormalized  version  of  e^n 
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A   normalized   backward   error  associated  with   the   prediction  of  y((m,n)-q) 
from  the  next  q-1  points  of  qLp.  is  defined  by 


rmqn-q  =   bA^(mn-q)yAl(k-l)y>-2) 


(6.9) 


As  with  the  forward  error  prediction  coefficients,  the  bAqA1(mn-q)  can  be 
interpreted  as  components  of  a  second  order  covariant  tensor.  The  components 
bAqV(mn-q)  equal  zero  for  the  range  of  indices  (A^Aj)  <  (mn-q)  or  (AX,A2)  >  (m,n). 
For  the  case  when  the  index  (Aj.A^  =  (mn-q)  the  component 

1 


bmqniq((m,n)-q)  = 


I  rq_1    I 


(6.10) 


where  r£~lq  is  an  unnormalized  version  of  F^q. 
1.     Normalized  Order  Update  Recursions 

The  following  two  theorems  are  stated  without  proof.  The  proofs  are 
identical  to  those  of  Theorems  5.1  and  5.2  of  the  previous  chapter  and  so  are  not 
repeated. 

a.     Theorem  6.1:  Normalized  Nonlinear  Levinson  Algorithm 

The  nonlinear  prediction  error  coefficients  can  be  updated  in  order 
using  the  following  recursions 
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b.     Theorem  6.2:  Normalized  Error  Order  Update  Algorithm 

The  nonlinear  prediction  errors  can  be  updated  in  order  according  to 
the  following  equation 
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(6.14) 


In  order  to  introduce  the  sample  y(k)  into  the  model  we  recognize 
the  backwards  errors  as  an  alternate  coordinate  system.  In  particular  the 
following  vector  is  defined. 
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Equations  (6.14)  correspond  to  a  model  structure  identical  with  that 
given  in  Figure  6.2.  In  this  diagram,  the  backwards  errors  given  in  equations 
(6.15)  correspond  to  those  that  are  shown  leaving  the  uppermost  stages  of  the 
filter. 

The  forward  error  in  predicting  y(k)  from  the  Y(k)  tensor  is  given  by 


(p-i) 


(6.16) 


ek-    -'    =  y(k)  -  K,V 

Because  the  backward  errors  are  uncorrected  it  is  a  straight  forward 
matter  to  find  the  (Fourier)  coefficients.  KA  .    They  are  given  by 
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The  m-th  component  is 
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Km  =  E|y(k)r(^o)-m|  (6.18a) 

=  E|ekmr^.o)-m|  (6.18b) 

=    |    ;  ekm|    i  Ek-r^J  (6.18c) 

=    I    !ekm|    i  p*  (6.18d) 

A  normalized  version  of  the  forward  error  defined  in  (6.16)  is  thus  given  by 
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I     I    ek     !     I  A'  =  0 

Recognizing  that  y(k)  is  a  zero  order  forward  error  we  can  write  a 
normalized  version  of  the  q-th  order  forward  error  as 

ekq  =  — —q ek«-'   -   p*  r(?ij.J  (6.20) 

i       ekq 

This  last  result  follows  if  (6.19)  is  iterated,  calculating  successively 
higher  order  errors  or  from  equation  (4.73)  where  the  equivalence  of  the  lattice 
model  and  Fourier  series  was  established. 

We  conclude  that  expressing  y(k)  as  a  stochastic  Fourier  series  (with 
the  backward  errors  as  a  basis)  amounts  to  nothing  more  than  the  addition  of  a 
supplementary  channel  to  the  existing  lattice  structure.  Equations  (6.11)  and 
(6.14)  hold  for  thi>  additional  channel  and  can  thus  be  used  to  update  the 
normalized  error  signals  and  model  parameters.  Figure  6.1  illustrates  a  quadratic 
nonlinear  lattice  filter. 

2.      Uniqueness  Of  Lattice  Parameters 
a.      Theorem  6.3 

The  lattice  parameterization  of  the  sequence  y(k)  is  unique.  That  is. 
the  lattice  parameters  given  by  (6.13)  are  unique. 


148 


b.     Proof  (Outline) 

The  truth  of  this  theorem  is  a  consequence  of  the  fact  that  the  lattice 
can  be  interpreted  as  a  Fourier  series  expansion  of  the  sequence  y(k)  in  terms  of 
the  orthonormal  backwards  errors  (see  equation   (6.17).)     The  uniqueness  of  the 
Fourier  series  is  well  known  [Ref.  47]. 
3.     Synthesis  Models 

A  model  analogous  to  that  illustrated  in  Figure  4.7  can  be  used  to 
regenerate  the  original  sequence.  This  requires  knowledge  of  the  forward  error 
residual  sequence.  The  other  inputs  (zero  order  forward  errors)  required  can  all 
be  obtained  from  past  estimates  of  the  sequence  or  from  initial  conditions. 

The  probability  distributions  of  the  output  forward  error  sequences  must 
be  known  if  a  synthesis  model  is  to  be  constructed  which  accurately  reproduces 
the  original  signal  statistics.  It  is  not  clear  that  any  general  statements  can  be 
made  about  the  nature  of  these  distributions.  It  has  been  shown  [Ref.  48:  p.  357] 
that  in  the  continuous  case  they  will  be  gaussian  and  of  the  same  variance  as  the 
original  additive  gaussian  noise  source.  In  the  same  reference  it  is  also  shown 
that  for  the  discrete  case  the  error  sequence  will  not  be  gaussian  although  it  will 
be  white.  The  inaccuracies  introduced  through  the  use  of  gaussian  noise  need  to 
be  investigated.  The  stability  of  these  lattice  models  is  also  left  as  an  open 
question. 

B.    SIMULATION  RESULTS 

In  order  to  verify  the  theory,  a  computer  simulation  was  performed.  Uniform 
white  noise  was  used  to  excite  the  system.  This  choice  provides  a  bounded  input 
so  that  the  response  of  the  system  used  for  the  simulation  could  be  guaranteed  to 
be  stable  for  a  suitable  choice  of  system  parameters.  A  sufficient  condition  for 
the  AR  model  to  be  stable,  given  a  driving  signal  whose  absolute  value  is 
bounded  on  (-a.a).  where  0  <  a  <   1.  is  given  by 


£       '       E      Hv     V    +  I  ai    <  l  (6.21) 

A,=l  AN=1 

This  condition  guarantees  that  the  output  never  exceeds  unity  and  thus  remains 
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Figure  6.1:  Quadratic  Nonlinear  Lattice  Filter 
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bounded.  This  condition  is  extremely  restrictive,  however,  it  does  provide  a 
simple  rule  for  building  models  which  can  be  used  for  simulation  purposes.  (The 
issue  of  the  inaccuracy  introduced  by  using  uniform  noise  to  drive  the  system  has 
been  avoided.) 

The  model  tested  was 


(6.22) 


Using   the   nonlinear  form   of  the   Levinson   algorithm,   equation    (6.11)    the 
following  model  parameters  were  estimated  for  two  repetitions  of  the  experiment, 
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VII.    CONCLUSIONS  AND  DISCUSSION 

A.    SUMMARY  OF  NEW  RESULTS 

In  this  dissertation  the  use  of  tensor  concepts  in  the  field  of  signal  processing 
was  investigated.  The  research  was  successful  in  a  number  of  areas,  extending 
known  results  and  introducing  some  new  ones.  In  particular,  tensors  were  used 
to  study  nonlinear  and  multidimensional  systems. 

The  nonlinear  modelling  problem  was  formulated  using  tensors.  By 
interchanging  the  roles  played  by  the  time  order  and  the  nonlinear  order  a  new 
form,  different  from  (although  equivalent  to)  the  traditional  Volterra  Series  was 
developed.  Using  this  new  form,  tensor  equivalents  of  the  discrete-time  Wiener- 
Hopf  and  Yule-Walker  equations  were  derived.  These  equations  can  be  solved  for 
the  unknown  system  parameters.  Conditions  for  the  solution  of  the  Wiener-Hopf 
equations,  without  requiring  matrix  inversion,  were  specified.  This  resulted  in  a 
computational  saving  of  Offp^l)3^"1')  operations,  where  N  is  the  largest  time 
delay  present  and  p  is  the  highest  exponent  present  in  the  system  model.  It  was 
further  shown  that  linear  adaptive  algorithms,  such  as  LMS  and  RLS.  can  be 
applied  to  solve  for  the  system  parameters. 

Existing  Lattice  filter  algorithms  were  reformulated  in  a  tensor  framework. 
It  was  shown  that  they  can  be  considered  to  be  equivalent  expansions  in 
alternate  coordinate  systems.  These  results  were  then  applied  to  the  solution  of 
the  2-D  autoregressive  modelling  problem.  Several  new  2-D  lattice  structures 
were  deduced.  These  models  are  not  efficient  in  the  sense  that  an  AR  model 
possessing  0(N'J)  parameters  would  require  0(N4)  parameters  when  recast  into  a 
lattice  form.  It  was  shown  that  with  proper  assumptions  of  shift-invariance  the 
complexity  of  the  lattice  models  can  be  reduced  to  0(N")  parameters. 

The  2-D  lattice  models  developed  were  then  used  to  deduce  a  nonlinear 
lattice  model.  This  model  differs  from  that  of  other  researchers  in  that  it  allows 
updates  to  be  performed  in  the  nonlinear  order  as  well  as  time  order. 
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Finally,  it  was  shown  that  these  lattice  models  are  amenable  to  systolic  array 
implementations. 

B.    FUTURE  DIRECTIONS  FOR  RESEARCH 

The  multidimensional  and  nonlinear  lattice  theories  are  by  no  means 
complete.  So  far,  several  new  lattice  models  have  been  developed.  It  now 
remains  to  investigate  the  properties  of  the  models. 

For  the  2-D  lattice  filters,  conditions  for  model  stability  have  yet  to  be 
established.  The  stability  of  multidimensional  systems  is  a  much  more  difficult 
problem  than  for  the  1-D  case  because  it  is  usually  impossible  to  factor 
multivariate  polynomials.  However,  necessary  and  sufficient  conditions  for  AR 
2-D  model  stability  have  been  established.  It  is  believed  that  these  can  be  used 
to  derive  lattice  stability  conditions. 

Another  important  topic  is  the  synthesis  of  2-D  transfer  functions  using 
lattice  structures.  Stated  differently,  the  problem  is  to  design  a  2-D  lattice 
parameter  filter  that  will  implement  a  given  2-D  frequency  response 
characteristic. 

For  the  nonlinear  lattice  structures  similar  issues  need  to  be  addressed.  In 
this  case,  however,  the  problems  are  more  complex  since  the  stability  and 
synthesis  questions  have  not  been  solved  for  the  nonlinear  AR  models. 

Before  the  synthesis  problem  for  the  lattice  implementation  can  be  solved 
several  more  basic  questions  need  to  be  answered.  The  first  is  the  definition  of  a 
useful  and  meaningful  specification  of  the  desired  filter  characteristic.  The 
nonlinear  AR  filter  affects  not  only  the  frequency  content  of  the  driving  signal 
but  also  its  probability  distribution.  The  filter  will  also  have  different  frequency 
characteristics  for  different  driving  signals.  All  of  these  effects  should  be  included 
in  the  specification. 

It  ha-  been  shown  that  for  the  continuous  case,  a  nonlinear  whitening  filter 
will  yield  a  signal  with  a  gaussian  probability  distribution  [Ref.  48:  p.  357j.  For 
the  discrete  case  it  has  been  shown  that  this  is  not  true.  The  inaccuracy 
introduced   by    approximating   the    input    to   the   synthesis   model    by    a    gaussian 
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introduced   by   approximating   the   input    to   the   synthesis   model  by   a   gaussian 
signal  is  an  important  question. 

Stability  of  the  nonlinear  lattices  is  also  an  important  question.  No  simple 
stability  tests  exist  for  the  nonlinear  AR  models.  Their  recursive  nature  makes 
stability  difficult  to  analize.  In  general,  stability  will  be  a  function  of  the  input 
as  well  as  the  system  which  significantly  complicates  the  problem. 
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APPENDIX  A:  ALTERNATE  PROOF  OF  THEOREM  4.2 

In  this  appendix  an  alternate  proof  of  the  generalized  error  order  update 
equations  (4.40)  is  provided.  The  proof  presented  here  relies  on  geometric 
arguments  which  are  possible  if  an  abstract  mathematical  framework  is  adopted. 
In  this  appendix,  a  random  process  will  be  considered  to  consist  of  a  time-series 
of  random  variables.  Each  of  these  random  variables  is  thought  of  as  a  vector  in 
an  infinite  dimensional  Hilbert  space.  For  a  rigorous  discussion  of  these  concepts 
see  for  example  [Ref.  51,52]. 

A.    DEFINITIONS  AND  FORMULATION 

A   discrete  random  process  Y  =  w(i):  0  <  i  <  K>  can  be  considered  to  span  a 

K+l  dimensional  subspace  Sk  +  1  of  an  infinite  dimensional  Hilbert  space  S 
(assuming  completeness  and  the  linear  independence  of  the  y(i)).  The  inner 
product  on  this  space  is  defined  to  be 

<u,v>  =  e{uv>  (A.l) 

where  u  and  v  are  elements  of  S.    This  inner  product  induces  a  norm 

i  I  u      =  <u.u>l; '  (A-2) 

The  error,  ekNi,.  in  predicting  the  element  y(k  +  l)  from  the  previous  N 
elements  of  Y  is  given  by 

e£,  =    £  hAN(k-l)y(A)  (A.3) 

where 

hJMk-1)  0    ■■■    (i       htNN.,        hkNN.2      •  hkN     1    0    •••    0:  (A. 4) 

A  normalized  version  of  the  forward  error  is  given  by 

N 

ekN_,  =  ^ (A. 5a) 


e 
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£>AN(k-l)y(A)  (A. 5b) 


where 


aAN(k+l)  =  ^  hAN(k+l)  (A.6) 

The    backwards    prediction    error,    rk^.N,    is    the    error    associated    with    the 
prediction  of  y(i)  given  the  next  N  elements  of  Y.    It  is  given  by, 


rk%  -    £  hAN(k-N)y(A)  (A. 7) 

where 

fhjftk-N)]  =  10   •  •  •   0   1    -hkNN+1    -hkN_N+2   •  ■  •    -hkN   o   ■  •  ■   0  (A. 8) 

A  normalized  version  of  r^.B  can  be  defined  as 


N 

r"kNN  =  ^r — r  (A. 9a) 

rk-  N 


-N  fk-N 


K 

Ni 


=    E  bAN(k-N)y(A)  (A.9b) 

A=0 


where 


hi"(k-N) 
bjf(k-N)  =        x\      '  (A. 10) 

Tk-N  I     ! 

By  orthogonal  we  mean  that  given  two  elements  of  the  space  S.  u  and  v.  say. 
then 

(0  for  u  ^  v  ,  .        . 

f  A. 11) 

u  lor  u  =  v 

This  will  be  indicated  symbolically  as 

u_v  (A.  12) 

The  expression 
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gM 


;a.i3) 


implies  that  u  is  orthogonal  to  all  elements  of  the  subspace  SM. 
The  symbol  '  ©"  is  used  to  indicate  direct  sum.    For  example 

S3  =  S1   ©  S2  (A. 14) 

means  that  the  space  spanned  by  S3  is  the  the  space  spanned  by  the  Cartesian 
product  of  the  underlying  sets  of  the  spaces  S1  and  S2  [Ref.    47:  p.  196]. 

Finally,  the  symbol  'v  '  will  be  used  to  mean  the  span  of. 

These  definitions  are  sufficient  to  state  and  prove  the  generalized  error  order 
update  theorem. 

B.    ERROR  ORDER  UPDATE  RECURSIONS 
1.     Theorem  4.2 

The   (N  +  l)   order  errors  can  be  calculated  from  the   N-  th  order  errors 
through  the  recursion 


-N  +  l 
ek^l 

—N  +  l 

rk-N 


\/l  -   (PN 


k    :)2 


i  „  k 

Pn  +  i  1 


-N 

ek  +  l 

-N 
rk-N 


'A.l, 


where 


k  ^— N     —  N      . 

PN-i-1    ~    <ek  +  l-rk-N> 


(A. 16] 


is  called  the  partial  correlation  coefficient  or  reflection  factor   [Ref.   17:  p. 


2.      Proof 


We  have  that 


fkNN  _  S=  *  |y(k-  N-  1)    y(k     N-2)  >(k) 


A.l' 


but 


rk    N        Skw+1   =    My(k-N)     y(k     Ni  1 


y(ki 


A.ll 


Therefore. 


1G2 


SkN^  =  v|FkNN   y(k-N-l)    -•■   y(k)|  (A. 19) 

which  can  be  written  as 

SkN+1  =  SkN   ©  v  L_\  (A.20) 

For  the  updated  forward  error,  eJVi1,  the  following  is  true 
e^l1  L  SkN  _  v  |FkN  \  (A.21) 

Since  e^,  j_  SkN   we   need   only  orthogonalize  it    (by   a   Gram-Schmidt   procedure) 
with  respect  to  v  ]rJlN>  to  obtain  e^t1-    Thus. 

•fit'  =  ekN+1  -  CrkN_N  (A.22) 

where 

N      — N 

<ek'_,,rk_N> 
C=  — _    k  N  (A. 23a) 

rk"-N 

-  P&+i       ekV    I  (A. 23b) 

Therefore. 

N 

-N+l  ek'l  r— N  k         — N  /    »      r>  4\ 

ek-l      =    ~ eK-l    "    />N-1   rk-N  (A. 24) 

ek-l 

Proof  of  The  fact  that 


N 
k-1 


*F-V  \/i  :  (pLu 


A.2I 


will  not  be  repeaTed  here  (see  Chapter  4.  Section  B). 

Similar  arguments  can  be  used  to  verify  The  backward  error  order  update 
equation  given  in  (A. 15). 


JG3 


APPENDIX  B:  FORTRAN  PROGRAM  LISTINGS 


c 

C        2DLAT  MAIN  PROGRAM 

C 

C       PROGRAM  TO  TEST  2D  LATTICE  SUBROUTINES 

C 

C        WRITTEN  29  APRIL  1985 

C 

c 

REAL*8  ALPHA(26,26),BETA(26,26),A(26,26),B(26,26),RHO(26,26) 

REAL*8  HZ(25).Y(256.256),R(26,26) 

REAL*4  GAIN(40,40) 
C 

C        DEFINE  THE  AR  FILTER  COEFFICIENTS 
C 

C        DATA  HZ/1.0,-  23, .12,  111, 21*0.0/ 
C        DATA  HZ  '1.0,-. 470, .007,.295,0.0,.015,. 055, .003,.022, 16*0.0/ 

DATA  HZ   1.0. .03.-  015,-.011,.033,-.47.. 195. .055.0.0,-. 085,0.0.0.0, 

*.  003.  .022.  -.0001,  .0067, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,-.  002, 0.0,  .0001/ 

c 

C        DEFINE  OTHER  PARAMETERS 

C 

C        MODEL  SIZE 

C 

N  =  5 

MN  =  N*N 
C 

C        ACTUAL  SYSTEM  SIZE 
C 

NA  =  5 

MNA  =  NA*NA 
C 

C        NUMBER  OF  POINTS  IN  THE  PLOT  IS  IP  X  IP. 
C 

IP  =  40 
C 

C        NUMBER  OF  DATA  POINTS  TO  ESTIMATE  CORRELATIONS  IS  IYS  X  IYS. 
C 

IYS    -   128 
C 

C        CALL  APPROPRIATE  SUBROUTINES 
C 

CALL  TXFCN(HZ.NA.GAIN\IP.0) 

CALL  PLOT8(GAIN.lP.IP. 'ORIGINAL  SPECTRUMS') 

CALL  GENRAT(HZ.NA.Y,IYS) 

CALL  CORLAT(Y,IYS,R,N) 

CALL  SCHURfRHO.R. ALPHA, BETA. MN) 

CALL  LEYSON(A.B.RHO.R.MN) 
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DO  10  I  =  1.  MN 
HZ(I)  =  A(1.I)/A(I,I) 
10      CONTINUE 

CALL  TXFCN(HZ,N.GAIN.IP.l) 

CALL  PLOT3(GAIN.IP.IP/4-TH  ORDER  LATTICE  APPROXIMATIONS- 
STOP 
END 
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C 

C        SUBROUTINE  GENRAT 

C 

C       THIS  SUBROUTINE  GENERATES  A  2D  DATA  FIELD  FROM  THE  AR 

C        FILTER  COEFFICIENTS.    A  WHITE  NOISE  INPUT  AND  WHITE  NOISE 

C        INITIAL  CONDITIONS  ARE  USED. 

C 

C        WRITTEN  29  APRIL  1985 

C 

C 

SUBROUTINE  GENRAT(HZ,N,Y,IYS) 
REAL*8  HZ(25),Y(256,256),ADD 
C 

C        FETCH  THE  RANDOM  NUMBER  GENERATOR  SEED 
C        UNDER  FORTHX.  STORE  SEED  IN  FILE  13. 
C        UNDER  DISSPLA  ENTER  SEED  EXPLICITLY. 
C 

C         READ(13.10)  IY 
CIO      FORMAT(IIO) 
C         REWIND  13 

IY  =  817346 
C 

C        GENERATE  THE  RANDOM  FIELD 
C 

MN  =  N*N 
MNM1  =  MN-1 
DO  40  1  =  l.IYS 
DO  30  J  -  l.IYS 

Y(I.J)  =  URAXD(IY)  -  .5 

KI  =  0 

KJ  -  0 

DO  20  K  -  1.MNM1 

IF  (MOD(K.N)  NE.O)  GOTO  15 
Kl  -  KI  -  1 
KJ       0 
GO  TO  16 

15  KJ    -   KJ  -    1 

16  IF  (((I-KI).LE.O)  OR.(fJ-KJ).LE.O))  GOTO  17 
ADD       ^  (I-KI.J-KJ) 

GO  TO  18 

17  ADD  -    URAND(IY)  -  .5 

18  Y(l.J)  =  Y(I.J)  -  HZ(K  -  1)*ADD 
20              CONTINUE 

30  CONTINUE 

40      CONTINUE 

(' 

('        STORF  THE  RANDOM  NUMBER  SEED 

C 

C         WRITE(  13.50)  IY 

C50      FORMAT(IIO) 

C 
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c 

RETURN 
END 
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C 

C  SUBROUTINE  CORLAT 

C 

C  THIS  SUBROUTINE  PRODUCES  A  CORRELATION  MATRIX  FROM 

C  A  2-D  DATA  FIELD  IN  AN  ORDER 

C  WHICH  IS  COMPATIBLE  WITH  SUBROUTINE  SCHUR 

C 

C  WRITTEN  30  APRIL  1985 

C 

P******************************************************  ***************** 

C 

SUBROUTINE  CORLAT  (Y,IYS,R,N) 
REAL*8  Y(256,256),R(26,26),SUM 
C 

C        DEFINE  CONSTANTS 
C 

MN  =  N*N 
IR  =  0 
C 

C        BEGIN  OUTER  LOOP 
C 

DO  100  MP1  =  l.N 
MO  =  MPl  -  1 
LLIM  =  2*MP1  -  1 
DO  90  L  =  1,LLIM 
L0  =  L  -  1 
II  =  MO 
Jl  =  L0/2 

IF  (MOD(L0.2).EQ.O)  GO  TO  10 
II  =  Jl 
Jl  =  M0 
10  JR  -  IR 

IR  =  IR  -  1 
KBOT  =  L 
DO  80  NPl  =  MP1.N 
NO  =  NPl  -  1 
KLIM  -  2* NPl  -  1 
DO  70  K  =  KBOT. KLIM 
K0   =  K-l 
I  =  NO 
J  =  K0/2 

IF  (MOD(K0.2).EQ.0)  GO  TO  20 
1    -  J 
J  =  NO 
20  IOFF  =  1-11 

JOFF  =  J- J] 
K1MAX  =  IYS 
K1MIN  =  IOFF  -    1 
IF  (IOFF.GT.0)  GO  TO  30 
KlMAX  =  IYS  +  IOFF 
KlMIN  =  1 
30  K2MAX  =  IYS 
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K2MIN  =  JOFF  -r  l 
IF  (JOFF.GT.O)  GO  TO  40 
K2MAX  =  IYS  -  JOFF 
K2MIN  =  1 
40  SUM  =  0.0 

JR  =  JR  +  1 
C  WRITE(6,44)IR,JR,Il,Jl)I,J,IOFF,JOFF 

C44  F0RMAT(8(2X,I4)) 

DO  60  Kl  =  K1MIN,K1MAX 
DO  50  K2  =  K2MIN,K2MAX 

SUM  =  SUM  +  Y(Kl,K2)*Y(Kl-IOFF.K2-JOFF) 
50  CONTINUE 

60  CONTINUE 

R(IR,JR)  =  SUM/(K1MAX-K1MIN+1)/(K2MAX-K2MIN-1) 
70  CONTINUE 

KBOT  =  1 
80  CONTINUE 

90  CONTINUE 

100     CONTINUE 
C 

C       FILL  IN  THE  SYMMETRIC  HALF  OF  CORRELATION  MATRIX 
C 

DO  120  I  -  2.MN 
1M1  =  1-1 
DO  110  J  =  1,IM1 
R(l.J)  =  R(J,I) 
110         CONTINUE 
120     CONTINUE 
C 
C 

RETURN 
END 
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C 

C  SUBROUTINE  TXFCN 

C 

C  GENERATES  A  2-D  FREQUENCY  RESPONSE  GIVEN  THE  TRANSFER 

C  FUNCTION  COEFFICIENTS 

C 

C  WRITTEN  BY  DR.  B.  MADAN 

C  MODIFIED  BY  P.J.  LENK    29  APRIL  1985 

C 

C  PARAMETER  KO  INDICATES  THE  ORDERING  OF  THE  COEFFICIENTS 

C  -  KO  =  0:  PARAMETERS  IN  ROW-MAJOR  ORDER 

C  -  KO  =  1:  PARAMETERS  IN  TWIDDLED  ORDER 

C 

p* *******************************************************  *************** 

C 

SUBROUTINE  TXFCN(HZ,N,GAIN,IP,KO) 

REAL*8  HZ(25),CONVRT(5,5) 
REAL*4  GAIN(40,40) 
COMPLEX  CSUM 
COMPLEX*8  ARG1.ARG2 
C 

C        DEFINE  CONSTANTS 
C 

PI  =  3.14159 
C 

C        DETERMINE  ORDER  AND  REORDER  IF  NECESSARY  THE  COEFFICIENTS 
C 

IF  (K0.EQ.0)  GO  TO  60 
JR  -  0 

DO  30  MPl  -   1,N 
MO  -  MPl  -  1 
LLIM  =  2*MP1  -  1 
DO  20  L  =  l.LLIM 
L0  =  L  -  1 
I  =  MO 
J  =  L0/2 

IF  (MOD(L0.2).EQ.O)  GO  TO  10 
I  =  J 
J  =  M0 
10  JR  --  JR  -    1 

CONVRT(]  -  l.J  +  1)  -  HZ(JR) 
20  CONTINUE 

30      CONTINUE 
( 

C        TRANSFER  THE  COEFFICIENTS  BACK  TO  HZ 
C 

JR  =  0 

DO  50  I  =   l.N 
DO  40  J  =  FN 
JR  -  JR  -  1 
HZ(JR)  =  CONVRT(FJ) 
WRITE(6.35)HZ(JR) 
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35  F0RMAT(2X.E12.5) 

40  CONTINUE 

50      CONTINUE 
C 

C       PROCEED  WITH  TRANSFER  FUNCTION  EVALUATION 
C 

60      DW  =  2.0*PI/FLOAT(IP-1) 
DO  100  IW1  -  1,IP 

Wl  =  DW  *  (IW1  -  1)-  PI 
DO  90  JW1  =  1,IP 
Al  =  Wl 

W2  =  DW  *  (JW1  -  1)-  PI 
IZ  =  0 

CSUM  =  CMPLX(0.0,0.0) 
DO  80  I  =  1,N 

ARG1  =  CMPLX(0  0,-Al) 
Al  =  Al  +  Wl 
A2  =  W2 
DO  70  J  =  l.N 
IZ  =  IZ  +  1 

ARG2  =  CMPLX(0.0,-A2) 
A2  =  A2  +  W2 
C  WRITE(6,77)I,J,ARG1,ARG2,1Z,HZ(IZ),CSUM 

C77  F0RMAT(2(2X!I3)!4(2X,E12.5),2X,I3,3(2X,E12.5)) 

CSUM  =  CSUM  +  CMPLX(SNGL(HZ(IZ)).0.0)*CEXP(ARGl-ARG2) 
70  •  CONTINUE 

80  CONTINUE 

GAIN(IWUWl)  -  l.O'CABS(CSUM) 
GAIN(IWi.JWi)  =  GAIN(IW1,JW1)*GAIN(IW1,JW1) 
C  WR1TE(6.78)IW1.JW1.CSUM.GAIN(IW1,JW1) 

C78  F0RMAT(2(2X.13).3(2X.E12.5)) 

90  CONTINUE 

C  \YRITE(6J8)IWl:(GAIN(IWl.I),I  =  l,IP) 

C18  F0RMAT(1X.I3,5(2X.E12.5)) 

100     CONTINUE 
RETURN" 
END 
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c 

C  SUBROUTINE  PLOTS 

C 

C  THIS  IS  A  ROUTINE  TO  USE  THE  DISSPLA  PACKAGE  TO 

C  DRAW  A  THREE  DIMENSIONAL  PLOT  OF  A  2-D  FILTER'S 

C  FREQUENCY  RESPONSE. 

C 

C  WRITTEN  BY  DR.  B.  MADAN 

C  MODIFIED  BY  P.J.  LENK    29  APRIL  1985 

C 

P******************************** *************************************** 

C 

SUBROUTINE  PLOT3(A,IM,JN,LABEL) 

DIMENSION  A(IM,JN) 

INTEGER  LABEL(IOO) 

CALL  TEK618 
C 

C        INITIALIZE  THE  PLOTTING  DEVICE 
C 

WR1TE(6.51) 

51  FORMATf  1  CALL  TEK618  OK') 
C  CALL  THE  DEVICE 

CALL  RESET("ALL") 
WR1TE(6.52) 

52  FORMATf  2  RESET  ALL  OK') 
C  CALL  SETUP  FOR  CUBE 

Al-FLOAT(IM) 
A2=FLOAT(JN) 
CALL  PAGE(11..9.5) 
\YRITE(6.53) 

53  FORMATC  3  CALL  PAGE  OK) 
C        CALL  PAGE(A1,A2) 

C        CALL  PHYSOR(0.,0.) 
CALL  AREA2D(7. 0.7.0) 
WRITE(6.54) 

54  FORMATf  4  CALL  AREA  2-D  OK") 
CALL  SIMPLX 

CALL  HEIGHT(.2) 

CALL  HEADIN(LABEL.lOO.l.l) 

CALL  HEJGHT(0.14) 

CALL  VIEW  (-10.0.-5  0.15.0) 

CALL  VOLM3D(12. 0.12.0.12.0) 
(  CALL  X  AXIS  LABELLING  ROUTINE 

CAIJ.  X3XAMEC  W2  8". 3) 
C  CALL  >    AXIS  LABELLING  ROUTINE 

CALL  Y3NAME(    W]  S'.SJ 
C  CALL  Z  AXIS  LABELLING  ROUTINE 

(ALL  Z3NAME("\2) 
C  CALL  THE  SURFACE  PLOT  ROUTINE 

C 

CALL  GRAF3D(-1. 0.0. 2. 1.0.- 1.0. 0.2. 1.0.0. 0.2. 0.16.) 
1500    CALL  SURMAT(A.l.IM.l.JN.O) 
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CALL  ENDPL  (1) 
CALL  DONEPL 
RETURN 
END 
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^******x****************xx***±*x»x*x*****xx:<:xvlxxx*x**Xx********x****X*xx 

c 

C         SUBROUTINE  SCHUR 

C 

C         CALCULATES  THE  REFLECTION  FACTORS  FROM  THE  CORRELATION  MATRIX 

C 

C         WRITTEN  29  APRIL  1985 

C 

P*** ***********************************************************  ********* 

C 

SUBROUTINE  SCHUR(RHO,R, ALPHA, BETA.N) 
REAL*8  RHO(26,26),R(26,26),ALPHA(26,26),BETA(26,26),RNORM,T 
C 

C         INITIALIZE  THE  ALPHA  AND  BETA  ARRAYS 
C 

DO  10  1  =  1,N 
DO  5  J  =  l.N 

ALPHA(LJ)  =  R(I,J)/DSQRT(R(I,I)) 
BETA(l.J)  =  ALPHA(IJ) 
RHO(I.J)  =  0.0 
5  CONTINUE 

C  WRITE(16.7)(ALPHA(1.J),J  =  1.N) 

C7  F0RMAT(5(2X:E12.5)) 

10      CONTINUE 
C 

C        BEGIN  CALCULATING  THE  REFLECTION  FACTORS 
C 

DO  50  J  -  2.N 
NJ1  =  N  -  J  +  1 
DO  40  I  =  l.NJl 
Jll  =  J  ■+■  I  -  1 
IP1  =  I  ~  1 

RHO(I.JIl)  =  ALPHA(LJI1)/BETA(IP1.JI1) 
RNORM  =  DSQRTfl.O-  RHO(I,JIl)*RHO(I,JIl)) 
DO  30  K  =  l.N 
T  -  ALPHA(I.K) 

ALPHA(l.K)  -  (ALPHA(I.K)-RHO(l.. II lpBETA(IPl.K)) /RNORM 
BETA(I.K)  -  (BETA(IPl.K)-RHO(I.Jll)*T)   RNORM 
30  CONTINUE 

40  CONTINUE 

('  WR1TE(16.42)J,((ALPHA(LK).K-1.N);I    l.N) 

('42  FORMAT)   2X.I3:4(2X.E12.5)) 

C  WR]TE(16.42)J.((BETA(I.K).K=1,N),I=1,N) 

50      CONTINUE 
C 
(' 

RETURN 
END 
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c 

C        SUBROUTINE  LEVINSON 
C 

C         GENERATES  THE  AUTOREGRESSIVE  MODEL  COEFFICIENTS  GIVEN  THE 

C         REFLECTION  FACTORS 

C 

C         WRITTEN  29  APRIL  1985 

C 

£-!**  ************************************************************  ********* 

c 

SUBROUTINE  LEVSON  (A,B,RHO,R,N) 
REAL*8  A(26,26),B(26!26),RHO(26,26),R(26,26),RNORM,T 
C 

C        INITIALIZE  THE  A  AND  B  MATRICES 
C 

DO  20  1  =  1,N 
DO  10  J  =  1,N 
A(I,J)  -  0.0 
B(I,J)  =  0.0 
10  CONTINUE 

A(I,I)  =  1.0/DSQRT(R(I,I)) 
B(LI)  =  A(LI) 
20      CONTINUE 
C 

C        CALCULATE  THE  AR  PARAMETERS  USING  LEVINSON'S  ALGORITHM 
C 

DO  60  J  =  2.N 
NJ1  =  N  -  J  +  1 
DO  50  I  =  l.NJl 
Ul  =  I  -t-  J  -  1 

RN'ORM  =  DSQRT(1.0-  RHO(I,IJl)*RHO(I,IJl)) 
DO  40  K  =  I.IJ1 
T  ==  A(l.K) 

A(I.K)  =  (A(l.K)  -  RHO(I.IJl)*B(I+l.K))'RNORM 
B(I.K)  =  (B(l-l.K)  -  RHO(I.IJl)*T)   RNORM 
40  CONTINUE 

50  CONTINUE 

C  WRITE(16.77)J.((A(1,K),K=1.N),I=1.N) 

C  WR1TE(16,77)J,((B(I,K).K=1.N)J-1.N) 

C77  F0RMAT(/2X,13.4(2X,E12.5)) 

60      CONTINUE 
DO  66  1  =  1,N 

RN'ORM  -  A (1.1)    A(l.l) 
C  WR1TE(6.65)RN0RM 

('  WRITE(17.65)RNORM 

C6f.  F0RMAT(2X,E12.5) 

m      CONTINUE 
C 

c 

RETURN 
END 
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c 

C  FUNCTION  URAND 

C 

C  TAKEN  FROM  "COMPUTER  METHODS  FOR  MATHEMATICAL  COMPUTATIONS"  BY 

C  G.E.  FORSYTHE.  M.A.  MALCOLM.  AND  C.B   MOLER 

C  PRENTICE-HALL,  ENGLEWOOD  CLIFFS,  NJ.,  1977,  P.  246. 

C 

j-i*********************************************************************** 

C 

REAL  FUNCTION  URAND (IY) 

INTEGER  IA,IC,ITWO,M2,M,MIC 
C         URAND  IS  A  UNIFORM  RANDOM  NUMBER  GENERATOR  BASED  ON  THEORY  AND 
C     SUGGESTIONS  GIVEN  IN  D.E.  KNUTH  (1969),  VOL.  2.    THE  INTEGER  IY 
C     SHOULD  BE  INITIALIZED  TO  AN  ARBITRARY  INTEGER  PRIOR  TO  THE  FIRST 
C     CALL  TO  URAND.    THE  CALLING  PROGRAM  SHOULD  NOT  ALTER  THE  VALUE  OF  IY 
C     BETWEEN  SUBSEQUENT  CALLS  TO  URAND.    VALUES  OF  URAND  WILL  BE  RETURNED 
C     IN  THE  INTERVAL  (0,1). 
C 

DOUBLE  PRECISION  HALFM 

REAL  S 

DOUBLE  PRECISION  DATAN,DSQRT 

DATA  M2/0/.ITWO/2/ 

IF(M2.NE.O)  GO  TO  20 
C 

C        IF  FIRST  ENTRY,  COMPUTE  MACHINE  INTEGER  WORD  LENGTH 
C 

M  =  l 
10  M2  =  M 

M=ITWO*M2 

IF(M.GT.M2)  GO  TO  10 

HALFM=M2 
C 

C        COMPUTE  MULTIPLIER  AND  INCREMENT  FOR  LINEAR  CONGRUENTIAL  METHOD 
C 

IA=8*IDINT(HALFM*DATAN(1.D0)/8.D0)^5 

IC=--2*IDINT(HALFM*(.5D0-DSQRT(3.D0)/6.D0))  +  ] 

MJC-(M2-1C)-M2 
C 

C        S  IS  THE  SCALE  FACTOR  FOR  CONVERTING  TO  FLOATING  POINT 
(' 

S=.5  HALFM 
C 

C        COMPUTE  NEXT  RANDOM  NUMBER 
C 

20  IY  =  IY*1A 
C 

('        THE  FOLLOWING  STATEMENT  IS  FOR  COMPUTERS  WHICH  DO  NOT  ALLOW 
C        INTEGER  OVERFLOW  ON  ADDITION 
C 

IF(IY.GT.MIC)  IY=(IY-M2)-M2 
C 

IY    IY+IC 
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c 

C       THE  FOLLOWING  IS  FOR  COMPUTERS  FOR  WHICH  THE  WORD  LENGTH 

C       FOR  ADDITION  IS  GREATER  THAN  FOR  MULTIPLICATION 

C 

IF(IY/2.GT.M2)IY=(IY-M2)-M2 
C 

C       THE  FOLLOWING  STATEMENT  IS  FOR  COMPUTERS  WHERE  INTEGER  OVERFLOW 
C       AFFECTS  THE  SIGN  BIT 
C 

IF(IY.LT.0)IY=(IY+M2)+M2 

URAND=FLOAT(IY)*S 

RETURN 

END 
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f1********************************************xxx»*****x:r**>************* 

c 

C        PROGRAM  NLMAIN 

C 

C        THIS  IS  THE  PROGRAM  TO  TEST  THE  NONLINEAR  LATTICE  MODEL 

C 

C        WRITTEN  7  MAY  1985 

C 

p********************************  ********************:,;***  *************** 

c 

REAL*8  ALPHA(26,26),BETA(26,26),A(26,26),B(26,26),RHO(26,26) 

REAL*8  HZ(5,5),Y(10000),R(26,26) 
C 

C        DEFINE  SYSTEM  PARAMETERS 
C 

DATA  HZ/. 2. .02,-. 1.0. 0,0.0, .05, 0.09,. 22,0.0,0.0,-. 03, .001,0.02,0.0 

*,0. 0,10*0.0/ 

c 

C        DEFINE  SYSTEM  CONSTANTS 
C 

DO  5  1  =  1,5 

K  -  6-  I 

WRITE(6,4)(HZ(K,J),J=1,5) 

4  F0RMAT(5(2X,E12.5)) 

5  CONTINUE 
IYS  =  1000 
NA  =  3 

MNA  =  NA  *  NA 

MNAPl  =  MNA  +  1 
C 

C        DEFINE  MODEL  CONSTANTS 
C 

N  =  3 

MN  =  N  *  N 

MNPl  =  MN  +  1 
C 

C        CALL  SUBROUTINES 
C 

CALL  NLGEN(Y.N.HZ.IYS) 

CALL  NLCLAT(Y.IYS.R.N) 

DO  30  I       l.MNPl 

WR1TE(16.20)(R(1,J),J  =  1,MNP1) 
20  FORMAT!  \5(2X.E12.5)) 

30      CONTINUE 

CALL  S(  HI  R(R  HO. H.  ALPHA. BETA. MNPl) 

DO  60  1  =  l.MNPl 

WRITE(16.20)(RHO(L  J),.I  =  l.MNPl) 
6(i       CONTINUE 

CALL  LEVSO\(A.B.RHO.R,MNPl) 

DO  70  I  =  l.MNPl 

\VRITE(16.20)(A(I.J),J=1,MNP1) 
70      CONTINUE 
C 
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STOP 
END 
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c 

C  SUBROUTINE  NLGEN 

C 

C  THIS  SUBROUTINE  GENERATES  AN  OUTPUT  SEQUENCE  FROM  THE  SYSTEM 

C  DESCRIBED  BY  THE  MODEL  PARAMETERS  CONTAINED  IN  H(,).    IT 

C  USES  WHITE  NOISE  UNIFORM  ON  (-.5, .5)  TO  EXCITE  THE  SYSTEM. 

C  THE  INITIAL  CONDITIONS  ARE  ALSO  DRAWN  FROM  THIS  DISTRIBUTION. 

C 

C  WRITTEN  MAY  7  1985 

C 

C 

SUBROUTINE  NLGEN(Y,N,H,IYS) 

REAL*8  Y(10000),H(5,5),DSEED 
C 

C        FETCH  THE  RANDOM  SEED 
C 

READ(13,10)  IY 
10      FORMAT(llO) 

REWIND  13 
C        DSEED  =  DFLOAT(IY) 
C 

C        SET  UP  THE  INITIAL  CONDITIONS 
C 

Y(l)  =  2.*(URAND(IY)  -  .5) 

Y(2)  =  2.*(URAND(IY)  -  .5) 
C        CALL  GGNML(DSEEDJYS.Y) 
C 

C        CALCULATE  THE  REMAINING  VALUES  OF  THE  SEQUENCE 
C 

DO  40  I  =  3.IYS 
C  Y(I)  =  2.*(URAND(IY)  -  .5) 

C  Y(I)  =    2*Y(I) 

Y(l)  =  URAND(IY) 
C  DO  30  J  =  l.N 

C  JMl  =  J  -  1 

C  DO  20  K   -   l.N 

C  KMl  =  K  -  1 

C  Y(I)       Y(l)  -  HIJ.KJ'COORDIYU-lj.JMlpCOORDtYfl^J.KMl) 

C20  CONTINUE 

C30  CONTINUE 

40      CONTINUE 
( 

('        FINISH 
C 
C         IY  -  DINT(DSEED) 

WRITE(13.50)IY 
50      FORMAT(IlO) 

REWIND  13 

RETURN 

END 
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C 

C  SUBROUTINE  NLCLAT 

C 

C  THIS  SUBROUTINE  PRODUCES  A  CORRELATION  MATRIX  FROM  NONLINEAR 

C  TIME  SEQUENCE  IN  AN  ORDER  WHICH  IS  COMPATIBLE  WITH  SUBROUTINE 

C  SCHUR. 

C 

C  WRITTEN  7  MAY  1985 

C 

C 

SUBROUTINE  NLCLAT  (Y,IYS,R,N) 
REAL*8  Y(10000),R(26,26),SUM,VEC(26) 
C 

C       DEFINE  CONSTANTS 
C 

MN  =  N*N 
MNP1  =  MN  +  1 
IYSM2  =  IYS  -  2 
FIYSM2  =  FLOAT(IYSM2) 
C 

C       INITIALIZE  R  MATRIX  TO  ZERO 
C 

DO  20  I  -  1.MNP1 
DO  10  J  =  l.MNPl 
R(I.J)  =  0.0 
10  CONTINUE 

20      CONTINUE 
C 

C       BEGIN  OUTER  LOOP 
C 

DO  80  I  =  3. IYS 
IR  =  1 

VEC(IR)  =  Y(I) 
DO  50  MPl  =  l.N 
M0  --  MPl  -  1 
LLIM  --  2*MPl  -  1 
DO  40  L  -    1  LLIM 
L0       L  -  1 
11        M0 
Jl  ~-  L0  2 

IF  (MOD(L0.2).EQ.0)  GO  TO  .",0 
11  -   Jl 
Jl  =  M0 
30  IR  -    IR  -  1 

VEC(IR)  =  COORD(Y(I-l).Il)*COORD(Y(I-2),Jl) 
40  CONTINUE 

50  CONTINUE 

C 

C       CALCULATE  THE  CORRELATIONS 
C 

DO  70  J    =    l.MNPl 
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DO  60  K  =  J,MNP1 

R(J.K)  =  R(J.K)  +  VEC(J)*VEC(K) 
C  WRITE(6.12)VEC(J),VEC(K).R(J.K) 

C12  F0RMAT(3(2X,E12.5)) 

60  CONTINUE 

70         CONTINUE 
80      CONTINUE 
C 

C       DIVIDE  BY  THE  NUMBER  OF  DATA  ELEMENTS  CONSIDERED 
C 

DO  100  J  =  1.MNP1 
DO  90  K  =  J.MNP1 

R(J,K)  =  R(J,K)/FIYSM2 
90  CONTINUE 

100     CONTINUE 
C 

C       FILL  IN  THE  SYMMETRIC  HALF  OF  CORRELATION  MATRIX 
C 

DO  120  I  =  2.MNP1 
IM1  =  1  -  1 
DO  110  J  =  1.IM1 
R(l.J)  =  R(J.I) 
110         CONTINUE 
120     CONTINUE 
C 
C 

RETURN 
END 
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c 

C  NONLINEAR  WIENER  MODELLING  PROGRAM 

C  THIS  USES  FUNCTIONS  CONTAINED  IN  ROUTINE  COORD 

C 

C  MODIFIED  1  JAN  85,  14  JAN  85 

C 

C 

DOUBLE  PRECISION  A(25,26),Z(25),X(15000),Y(15000),VAR,XA,YA,YHAT 

DOUBLE  PRECISION  ZZ(25),XMl,SEED,STORE(25) 
C 

C        INPUT  SIMULATION  PARAMETERS 
C 

READ(13,77)IY 

REWIND  13 

WRITE{6,40) 

40  FORMAT(2X, 'MAGNITUDE  OF  NOISE') 
READ(5.41)VAR 

41  F0RMAT(F12.5) 
WRITE(6,42) 

42  FORMAT(2X, 'NUMBER  OF  POINTS') 
READ(5,43)N 

43  F0RMAT(I5) 
WRITE(6.44) 

44  FORMAT(2X.'MAXIMUM  POWER  OF  X  ') 
READ(5.45)LW 

45  FORMAT(Il) 
WRITE(12.46)N 
WRITE(6.46)N 

46  FORMAT(2X.THE  NUMBER  OF  POINTS  WAS  ',15) 
WRITE(12.47)LW 

WRITE(6.47)LW 

47  FORMAT(2X.'THE  MAXIMUM  POWER  OF  X  IS  ',11) 
WRITE(6  48)VAR.VAR 

WR1TE(12  48)YAR.YAR 

48  FORM AT(2X. THE  NOISE  IS  UNIFORM  FROM  -'.F9.5."  TO  +',F9.5) 
L\V  =  LW  -  1 

VAR  =  VAR*2.0 

DO  132  I  =   1.25 
STORE(I)  -  0.0 
132     CONTINUE 

X(l)  -  VAR*(URAND(IY)  -  .5) 

Y(l)  =  1.0 

DO  2  1-2.X 

X(I)  -  YAR*(URAND(1Y)  -  .5) 
Y(I)  =  UNKNOW(X(I),X(I-l)) 
2       CONTINUE 

CALL  NCRLAT(X,Y.A.LSQW,LW,N) 
C         WRITE(6,134) 
C         WR1TE(12.134) 
C134     FORMAT(/2X,'NO  INVERSION  SOLUTION) 

DO  142  I  -  1.  LSQW 


183 


ZZ(J)  =  A(J,LSQW  +  1)/A(J.J) 
STORE(J)  =  STORE(J)  +  ZZ(J)/25.0 
142     CONTINUE 
C         DO  133  J  =  l.LW 
C  IMIN  =  (J-1)*LW  f  1 

C  IMAX  =  J*LW 

C  WRITE(6,173)  (ZZ(I),I  =  IMIN, IMAX) 

C  WRITE(12,173)  (ZZ(I),I  =  IMIN, IMAX) 

C173         F0RMAT(5(2X,E12.5)) 
C133     CONTINUE 
C         CALL  SOLVE(A,Z,LSQW) 
C         WRITE(6,135) 
C         WRITE(12,135) 

C135     F0RMAT(/2X,'USING  FULL  MATRIX  INVERSION') 
C         DO  27  J  =  1,LW 
C  IMIN  =  (J-l)'LW  +  1 

C  IMAX  =  J*LW 

C  WRITE(6.11)  (Z(I),I  =  IMINJMAX) 

C  WRITE(12,11)  (Z(I).I  =  IMINJMAX) 

11  F0RMAT(5(2X,E12.5)) 

C27      CONTINUE 
131     CONTINUE 
WRITE(6.201) 
WRITE(12,201) 
201     FORMAT(  '2X.NO  INVERSION  SOLUTION  AVERAGED  OVER  25  RUNS' 
DO  200  J  =  1,LW 

IMIN  =  (J  -  lj'LW  +  1 
IMAX  ='  J  *  LW 

WRITE(6.11)  (STORE(I).I  =  IMINJMAX) 
WRITE(12.11)  (STORE(I).I  =  IMIN. IMAX) 
200     CONTINUE 
\VRITE(6,22) 
WRITE(  12,22) 

XM1  =  VAR*(URAND(IY)  -  .5) 
DO  73  I  =  1.10 

XA  =  VAR*(URAND(IY)  -  .5) 
YA  =  UNKNOW(XA.XMl) 
YHAT  =  0.0 
YYHAT  -  0.0 
DO  14  J       l.LW 
JM1  -  J  -  1 

X2  =  COORD(XMl.JMl) 
DO  15  K  =  l.LW 
KM]       K  -  1 
X!   =  COORD(XA.KMl) 
C  WRITE(6.52)Xl.X2 

C52  F0RMAT(2(2X.E12.5)) 

YHAT  =   YHAT   -   Xl  *  X2  *  Z((J-1)*LW    •   K) 
WHAT   -   YYHAT   -  Xl  *  X2  *  ZZ((J-1)*LW   -   K) 
15  CONTINUE 

14  CONTINUE 

ERROR  =  (YA  -  YHAT)   YA 
ZERROR       (YA  -  YYHAT)   YA 
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XM1  =  XA 

WRITE(6,17)  YA,YHAT,YYHAT, ERROR. ZERROR 
17  F0RMAT(5(2X.E12.5)) 

22  FORMAT(/6X.'YAM2X,'YHAT'10X.!YYHAT'.9X.'ERROR'19X.'ZERROR') 

WRITE(12,17)  YA.YHAT.YYHAT. ERROR, ZERROR 
73      CONTINUE 

WRITE(13,77)IY 
77      FORMAT(IlO) 
STOP 
END 
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C 

C       SUBROUTINE  UNKNOW 

C 

C        THIS  ROUTINE  DEFINES  THE  UNKNOWN  SYSTEM 

C 

C 

DOUBLE  PRECISION  FUNCTION  UNKNOW(X,Xl) 
DOUBLE  PRECISION  H(25),YHAT,XTl,XT2 
LW  =  5 


H(l)  =    2 

H(2)  =  -.4 

H(3)  =  .03 

H(4)  =  -.7 

H(5)  =  0.0 

H(6)  =  .5 

H(7)  =  .35 

H(8)  =  .11 

H(9)  =  .9 

H(10)  =  0.0 

H(ll)  =  .01 

H(12)  =  13 

H(13)  =  -.33 

H(14)  =  .7 

H(15)  =  0.0 

H(16)  =  .43 

H(17)  =  .81 

H(18)  -  -.05 

H(19)  =    4 

H(20)  =  0.0 

H(21)  =  0.0 

H(22)  =  0.0 

H(23)  =  0.0 

H(24)  =  0.0 

H(25)  =  0.0 

VHAT  =  0.0 

DO  14  J  -   l.LW 

JM1  =  J  -  1 

XT2  =  COORD(XUJMl) 

DO  15  K  -  1,1.  W 

KM!  -  K  -  1 

XTl  -    COORD(X.KMl) 

('                 WRITE(6.52)XTl.XT2 

C52              F0RMAT(2(2X.E12  ".)) 

YHAT       YHAT    •   XTl  "  XT2  ' 

HfJYirLW 

15          CONTINUE 

14       CONTINUE 

UNKNOW        YHAT 

RETURN 

END 
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p**********************************************x*****************x****** 

c 

C       FUNCTION  COORD 

C 

C       GENERATES  OUTPUT  OF  THE  FUNCTIONS  BEING  USED  AS  COORDINATES 

C 

C       CREATED  23  AUG  84 

C 

p*** *****************************  *************************************** 

C 

DOUBLE  PRECISION  FUNCTION  COORD (X,I) 
C 

C  USE  LEGENDRE  POLYNOMIALS 
C 

C  IF  (I.NE.O)  GO  TO  1 

C  Y  =  1.0 

C  GO  TO  30 

CI  IF  (I.NE.l)  GO  TO  2 

C  Y  =  1.732051*X 

C  GO  TO  30 

C2  IF  (I.NE.2)  GO  TO  3 

C  Y  =  3.354102*(X*X  -  173.) 

C  GO  TO  30 

C3  Y=  6.61438*(X*X*X-  3./5.*X) 
C 

C  USE  SIMPLE  POWER  SERIES  TYPE  POLINOMIALS 
C 

Y  =  1.0 

IF  (I.EQ.0)  GO  TO  30 

Y  =  X**I 

30      COORD  =  DBLE(Y) 
RETURN 
END 


187 


C 

C        SUBROUTINE  SOLVE 

C 

C        SUBROUTINE  TO  SOLVE  A  SYSTEM  OF  LINEAR  EQUATIONS 

C        USES  GAUSSIAN  ELIMINATION  WITH  PARTIAL  PIVOTING 

C 


SUBROUTINE  SOLVE(A,Z,K) 

DOUBLE  PRECISION  A(25,26),Z(25),Y,TEMP,LARGE 


C 


KMl  =  K  -  1 

KP1  =  K  +  1 

c 

c 

DO  10  L  =  1,KM1 

LP1  =  L  +  1 
DO  11  I  =  LP1,K 

LARGE  =  DABS(A(L,L)) 

LROW  =  L 

DO  12  M  =  LP1.K 

IF  (DABS(A(M.Lj).LE. LARGE)  GO  TO  12 
LARGE  =  DABS(A(M,L)) 
LROW  =  M 

12  CONTINUE 
C 

IF  (L.EQ.LROW)  GO  TO  13 
DO  14  M  =  L,KP1 
TEMP  -  A(L.M) 
A(L.M)  =  A(LROW.M) 
•      A(LROW.M)  =  TEMP 

14  CONTINUE 
C 

13  Y  -  A(I.L)    A(L.L) 
C 

DO  15  M  -   l.KPl 

A(I.M)  =  A(I,M)  -  A(L:M)*Y 

15  CONTINUE 
11           CONTINUE 
10      CONTINUE 

C 

c 

Z(K)       A(K.KPl)   A(K.K) 
DO  16  L  -    l.KMl 
I    =  K  -  L 
Z(l)  -    A(I.KPl) 
KMl   =  K  -  I 
DO  17  \1       I  .KMl 
J  -  K  -  M  -  1 
Z(I)  =  Z(I)  -  A(1.J)*Z(J) 
17  CONTINUE 

Z(I)  =  Z(I)   A(I. I) 
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16      CONTINUE 

C 

C 

RETURN 

END 
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prr******************************************************************** 

c 

C  NONLINEAR  CORRELATION  MATRIX  CALCULATOR 

C 

C  CREATED  23  AUG  84 

C  USES  FUNCTION  UNKNOW  TO  DETERMINE  FUNCTIONS  FOR  EXPANSION. 

C  X  -  INPUT  VECTOR 

C  Y  -  OUTPUT  VECTOR 

C  PHI  -  CORRELATION  MATRIX 

C  -  LAST  COLUMN  CONTAINS  INPUT/OUTPUT  CROSS-CORRELATION 

C 

P*** *****************************************************  ************* 

C 


c 
c 


SUBROUTINE  NCRLAT(X,Y,PHI.LSQW,LW,N) 
DOUBLE  PRECISION  PHI(25,26),X(1).Y(1),TX(5,5) 


NMl  =  N  -  1 
LSQW  -  LW  *  LW 
LSQWPl  =  LSQW  +  1 


c 

c 

DO  40  1=1. LSQW 

DO  39  J  =  l, LSQWPl 

PHI(I.J)  =  0.0 

39 

CONTINUE 

40 

CONTINUE 

C 

c 

DO  2  M  =  2.N 

MMl  =  M-1 

DO  3  1=1. LW 

IMl  =  I  -  1 

X2  =  COORD(X(MMl),IMl) 

DO  4  J=1,LW 

JMl  -  J  -  1 

XI  =  COORD(X(M),JMl) 

c 

WRITE(6.81)  XI. X2 

C81 

F0RMAT(2(2X.E12.5)) 

TX(l.J)  =  XI  *  X2 

C 

VVRITE(6.80)  X(K),X(MM1),LJ,TX(I,J) 

C80 

FORM  AT(2X.2(2X.El2.5).2(2X.12).El  2.5 

4 

CONTINUE 

.1 

CONTINUE 

DO  5  1-  l.LW 
DO  6  II     l.LW 

K   =  (I-1)*LW  -   II 

PH1(K.  LSQWPl)  =  PHI(K.  LSQWPl)  +  Y(M)*TX(1,II] 
DO  7  J  =  1.LW 
DO  8  JJ=1,LW 

KK  -  (J-l)TW  -  .].] 
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PHI(K,KK)  =  PHI(K.KK)  +  TX(1,II)*TX(J.JJ) 

8  CONTINUE 

7  CONTINUE 

6  CONTINUE 

5  CONTINUE 

2       CONTINUE 

C 

DO  31  I=1,LSQW 
DO  32  J=1,LSQWP1 

PHI(I,J)  =  PHI(I,J)/FLOAT(NMl) 

32  CONTINUE 

31      CONTINUE 

C 

C 

C         D0  17I=1,LSQW 

C  WRITE(11,16)  (PHI(I,J),  J  =  l.LSQW) 

C16  F0RMAT(/,9(2X,E12.5)) 

C17      CONTINUE 

WRITE(11.19)  (PHI(I,LSQWP1),I  =  l.LSQW) 

19      F0RMAT(/,',9(2X,E12.5)) 

C 

C 

RETURN- 
END 
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C 

C  NONLINEAR  WIENER  MODELLING  PROGRAM 

C  THIS  USES  THE  LMS  ADAPTIVE  ALGORITHM  14  FEB  84  (VALENTINE'S) 

C 

C  UNKNOWN  SYSTEM  DEFINED  IN  FUNCTION  UNKNOW 

C 

*-,******************************** *************************************** 

C 

DOUBLE  PRECISION  TX(5,5),H(5,5),VAR,X,XMl,Y,YHAT,SEED 

READ(13,77)IY 

REWIND  13 

WRITE(6,40) 

40  FORMAT(2X,'MAGNITUDE  OF  NOISE') 
READ(5,41)VAR 

41  F0RMAT(F12.5) 
WRITE(6,42) 

4?      FORMAT(2X,'NUMBER  OF  POINTS') 
READ(5,43)N 

43  F0RMAT(I5) 
WRITE(6,44) 

44  FORMAT(2X, 'MAXIMUM  POWER  OF  X  ') 
READ(5,45)LW 

45  FORMAT(Il) 
WRITE(6.7) 

7  FORMAT(2X.  CONVERGENCE  FACTOR') 
READ(6.8)U 

8  FORMAT(Fl25) 
WRITE(12,46)N 
WRITE(6.46)N 

46  FORMAT(2X.THE  NUMBER  OF  POINTS  WAS  ',15) 
WRITE(12,47)LW 

WRITE(6.47)LW 

47  FORMAT(2Xt'THE  MAXIMUM  POWER  OF  X  IS  \Il) 
WRITE(6.48)VAR.VAR 
WR1TE(12.48)VAR.VAR 

48  FORMAT(2X.THE  NOISE  IS  UNIFORM  FROM  -\F9.5.'  TO  V,F9.5) 
WRITE(6.49)U 

WRITE(12.49)U 

49  FORM AT(2X. THE  CONVERGENCE  FACTOR  WAS  \E12.5) 
L\\    -   l.W    -    1 

VAR  =   YAR*2  0 
XMl  -  VAR*(URAND(1Y)  -  .5) 
(' 

(       IMTI  \LIZE  THE  11  TENSOR 
(' 

DO  10  I       l.LW 
DO  9  J       l.LW 
H(1J)   -   0.0 

9  CONTINUE 

10  CONTINUE 
C 

C     BEGIN  THE  ITERATION  LOOP 


192 


c 

DO  2  K=2,N 

X  =  VAR*(URAND(IY)  -  .5) 
Y  =  UNKNOW(X.XMl) 
C  WRITE(6.97)I.X.Y 

C97  F0RMAT(2X.I5.2X,E12.5,2X.E12.5) 

YHAT  =  0.0 
C 

C      CALCULATE  THE  MEASUREMENT  TENSOR  AND  THE  OUPUT  ESTIMATE 
C 

DO  14  I  =  1,LW 
IM1  =  I  -  1 

X2  =  COORD(XMl,IMl) 
DO  15  J  =  1,LW 
JM1  =  J  -  1 
XI  =  COORD(X,JMl) 
C  WRITE(6.52)X1,X2 

C52  F0RMAT(2(2X,E12.5)) 

TX(I.J)  =  Xl*X2 

YHAT  =  YHAT  +  TX(I,J)  *  H(I,J) 
15  CONTINUE 

14  CONTINUE 

C 

C      CALCULATE  THE  NEW  VALUE  OF  THE  TENSOR 
C 

ERROR  =  Y  -  YHAT 
DO  5  I  =  l.LW 
DO  4  J  -  l.LW 

H(LJ)  -  H(IJ)  +  2.0*U*ERROR*TX(I,J) 

4  CONTINUE 

5  CONTINUE 
XM1  =  X 

C 

C   PRINT  THE  NEW  VALUE  OF  THE  TENSOR 

C 

KM1  =  K  -  1 

IF  (KM1.NE.(KM1  100)*100)GO  TO  2 

WRITE(6.13)KMl 

WRITE(12.13)KM1 
13  FORMAT(   2X.  ITERATION  NUMBER  M5) 

WRITE(G  135) 

WRITEf  12.135) 
135         FORMATC'X   THE  NEW  VALUE  OF  THE  H  TENSOR  IS) 

DO  27  1        l.LW 

WRITE(C.II)  |H(I. .])..}  -    l.LW) 
WR1TE(12.11)  (H(l.J).J       l.LW  ) 
11  F0RMAT(5(2X.E12.5)) 

27  CONTINUE 

WRITE(6.22) 

WRITE(  12.22) 

WRITE(6.17)  Y. YHAT. ERROR 
17  F0RMAT(5(2X.E12.5)) 

22  FORMAT(6X.:Y".12X.YHAT'10X.  ERROR) 
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WRITE(12,17)  Y.YHAT.ERROR 
2        CONTINUE 

WRITE(13,77)IY 
77       FORMAT(IIO) 

STOP 

END 
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